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1.  Introduction 

This  final  report  contains  a  summary  of  work  accomplished  on  O.  N.  R.  Contract  N00014- 
86-K-0370,  High-Resolution  Radar-Imaging ,  for  the  period  from  15  January  1989  to  14  January 
1990.  Also  included  is  a  description  of  research  in  progress  that  will  be  phased  out  as  funding  for 
this  project  has  been  terminated. 

/The  goal  of  this  project  is  to  formulate  and  investigate  new  approaches  for  forming  images 
of  radar  targets  from  spotlight-mode,  delay-doppler  measurements.  These  measurements  can  be 
acquired  with  a  high-resolution  radar-imaging  system  operating  with  an  optical-  or  radio-frequency 
carrier.  Work  in  this  reporting  period  has  concentrated  on  our  estimation-theory  approach  to 
forming  high  resolution  images.  This  approach  accounts  for  measurement  noise  and  for  the  statistical 
properties  of  radar-backscatter  data.  /V  - 

2.  Summary  of  Work  Accomplished 

2.1  Conventional  Approach  to  Imaging 

A  computer  implementation  of  the  conventional  method  for  forming  radar  images  via  the 
two-dimensional  Fourier-transform  has  been  implemented  by  Mr.  D.  Porter,  who  is  an 
undergraduate  student  in  Electrical  Engineering.  This  is  used  to  compare  and  evaluate  images 
produced  conventionally  with  those  produced  using  the  new  methods  we  are  studying.  There 
are  two  modules  in  the  program. 

The  first  module  produces  simulated  radar  back-scatter  data.  The  simulation  calculates 
ideal  samples  of  a  received  signal  when  a  stepped-frequency  waveform  is  reflected  off  of  a 
target  having  a  specified  scattering  function  and  the  received  signal  is  mixed  with  the  transmitted 
signal  and  then  sampled  in  quadrature  to  produce  a  complex-valued  sample.  The  input  of  this 
simulation  includes  the  scattering  function  of  the  target,  the  desired  resolution,  and  the  base 
frequency  of  the  transmitted  waveform.  Complex  samples  of  the  received  signal  form  the 
output  of  the  program.  The  computational  complexity  is  0(N4),  where  N  is  the  number  of 
resolution  cells  in  one  target  dimension.  While  coding  efficiencies  have  been  exploited  to  a 
high  degree,  it  may  still  be  desirable  to  effect  a  parallel  implementation  of  the  calculations  for 
producing  results  more  quickly  for  high  resolution  images  when  N  is  large. 
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The  processing  of  simulated  and  real-data  samples  is  performed  by  the  second  module, 
which  applies  a  slightly  modified  two-dimensional  DFT  to  the  samples.  The  squared  magnitude 
of  the  resulting  transform  is  then  produced  as  the  target’s  image.  It  is  assumed  that  scatterers 
do  not  migrate  out  of  their  resolution  cells  over  the  time  of  data  collection,  which  is  an 
approximation.  When  this  approximation  is  not  well  met,  the  resulting  image  is  distorted  due 
to  range  walk.  A  number  of  ways  are  described  in  the  literature  to  deal  with  this  problem. 
First  of  all,  the  transmitted  waveform’s  base  frequency  can  be  chosen  sufficiently  high  so  that 
the  desired  cross  range  resolution  can  be  achieved  without  requiring  a  wide  variation  in  view 
angle  over  the  data  collection  interval.  Since  only  a  narrow  angle  is  required,  scatterers  will 
remain  within  the  confines  of  their  original  resolution  cell.  Secondly,  the  effect  of  range  walk 
can  be  seen  as  equivalent  to  collecting  polar-formatted  data  in  the  spatial  frequency  domain, 
yet  transforming  it  as  if  it  were  in  a  rectangular  format.  Interpolation  using,  for  example, 
cubic  splines  in  order  to  reposition  the  data  samples  onto  a  rectangular  grid  in  the  spatial 
frequency  domain  is  a  method  of  "focusing"  the  image,  thereby  reducing  the  error  due  to  range 
walk.  This  focusing  method,  which  we  are  using,  is  described  by  D.  Wehner  [High  Resolution 
Radar ,  Artech  House,  pp.  311-317,  1987], 

In  addition  to  the  two  main  modules  of  the  program,  additional  utility  modules  have  been 
developed  for  the  generation  of  scattering  function  files,  the  display  of  images  on  MASSCOMP 
and  SUN  workstations,  and  the  conversion  of  data  between  the  different  formats  used  on  the 
different  computers  used  in  our  study. 

2.2  Estimation-Theorv  Approach  to  Imaging 

Significant  progress  has  been  made  on  the  method  we  are  developing  for  producing  images 
of  low  visibility  targets  modeled  as  a  diffuse  scattering  object.  The  problem  of  estimating  the 
scattering  function  of  a  diffuse  target  is  ill  posed,  with  the  result  that  estimates  are  unstable 
having  a  the  rough  appearance  of  an  object  with  strong  speckle  noise.  Thus,  regularization  is 
required  in  order  to  stabilize  estimates  of  the  scattering  function,  which  is  a  two-dimensional 
power-density  spectrum.  Of  particular  importance  in  our  work  over  the  reporting  period  is 
the  identification  of  a  method  for  regularizing  estimates  via  the  method  of  sieves  introduced 
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by  U.  Grenander  ( Abstract  Inference,  Wiley,  1981).  We  expect  that  this  will  be  an  important 
development  not  just  for  radar  imaging  but  also  for  any  problem  where  a  power-density 
spectrum  must  be  estimated  from  noisy  data.  The  asymptotic  properties  of  this  method  of 
regularization  have  also  been  instrumental  in  the  identification  of  a  computational  method  for 
producing  target  images  practically.  The  inversion  of  large  matrices  is  not  required  under  some 
conditions  that  are  met  is  practice,  which  removes  a  major  impediment  previously  existing  with 
our  method. 

The  following  paragraphs  contain  the  abstracts  of  publications  describing  our  results. 
Details  are  given  in  the  appendices,  which  contain  reprints  of  the  publications. 

2.2.1  Abstract  of:  P.  Snvder.  J,  O’Sullivan,  and  M.  Miller.  "The  Use  of  Maximum  Like¬ 
lihood  Estimation  for  Forming  Images  of  Diffuse  Radar  Targets  from  Delav-Doppler 
Data."  IEEE  Trans,  on  Information  Theory.  Vol.  35.  dp.  536-548.  1989:  see  Appendix  6.1 
for  a  reprint. 

This  publication  gives  the  model  and  fundamental  estimation  equations  for  the  method 
we  are  developing.  The  abstract  is: 

"A  new  approach  to  high  resolution  radar  imaging  is  presented.  The  starting 
point  is  a  model  of  the  radar  echo  signal  based  on  the  physics  governing  radar 
reflections.  This  model  has  been  used  several  times  in  the  past  for  describing  radar 
targets  that  are  rough  compared  to  the  wavelength  of  the  transmitted  radiation. 
Without  specif ying  precisely  what  the  transmitted  signal  is,  a  general  estimation-based 
procedure  is  derived  for  obtaining  images.  After  discretizing  the  model,  the  radar 
imaging  problem  reduces  to  the  task  of  estimating  discretized  second-order  statistics 
of  the  reflectance  process  of  the  target.  Maximum  likelihood  estimates  of  these 
statistics  are  obtained  as  the  limit  point  of  an  expectation-maximization  algorithm." 
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2.2.2  Abstract  of:  P.  Moulin.  D.  Snvder.  and  J.  O’Sullivan.  "Maximum-Likelihood  Spec¬ 


trum  Estimation  of  Periodic  Processes  from  Noisv  Data."  Proc.  1989  CISS  Conference. 
Johns  Hopkins  University.  Baltimore.  MD.  March  1989:  see  Appendix  6.2  for  a  reprint. 

"We  have  developed  a  new  approach  to  maximum-likelihood  spectrum  estimation 
of  wide-sense  stationary  processes  from  noisy  data.  A  statistical  model  for  the  data 
is  defined.  The  process  whose  spectrum  is  sought  is  wide-sense  stationary,  periodic 
and  Gaussian,  and  its  observations  are  corrupted  by  an  additive  white  noise  [and  a 
linear  transformation],  [For  our  radar-imaging  problem,  the  Gaussian  process  models 
radar  back-scatter  data  from  a  diffuse  target,  the  spectrum  is  the  target’s  scattering 
function,  the  data  are  corrupted  by  additive  noise,  and  the  linear  transformation 
depends  on  the  transmitted  signal.]  A  maximum-likelihood  formulation  of  this 
problem  has  been  derived,  and  the  equations  are  solved  numerically  via  the 
expectation-maximization  algorithm.  This  approach  presents  several  attractive 
features,  an  important  one  being  that  the  noise  corrupting  the  observations  is  now 
taken  into  account. 

We  present  some  recent  developments  for  this  problem.  The  statistical  per¬ 
formance  of  the  new  maximum-likelihood  spectrum  estimator  is  studied  both 
theoretically  and  numerically.  Comparison  with  traditional  estimators,  such  as  the 
periodogram,  highlight  several  strong  points  of  the  method.  We  also  identify  certain 
limitations,  namely  the  instability  of  estimates  for  high  noise  levels,  [which  is  due 
to  the  ill-posed  nature  of  the  spectrum  estimation  problem].  These  limitations  can 
be  alleviated  if  a  priori  information  about  the  signal  is  available.  Two  such  problems 
are  discussed  [in  Appendix  6.2]  in  which  the  information  at  hand  has  the  form  of 
a  constraint  on  the  input  signal-to-noise  ratio. 

We  show  [in  Appendix  6.2]  how  such  information  can  be  incorporated  in  the 
maximum-likelihood  estimation  procedure.  First  we  assume  the  signal  power  to  be 
known.  Theoretical  issues  of  existence  and  uniqueness  of  the  solution  are  discussed. 
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We  proceed  with  a  problem  in  which  the  information  is  less  complete,  when  only 
an  upper  bound  and/or  a  lower  bound  on  the  signal  power  are  available.  The 
statistical  performance  of  both  constrained  estimators  is  quantitatively  studied.” 

2.2.3  Abstract  of:  J.  A.  O’Sullivan.  D.  L.  Snyder,  and  P.  Moulin:  "The  Role  of  Spectrum 
Estimation  in  Forming  High-Resolution  Radar  Images".  Proc.  ICASSP  1989.  Glasgow. 
U.K..  Mav  1989:  see  Appendix  6.3  for  a  reprint. 

"We  have  developed  a  new  approach  to  forming  high-resolution  images  of  radar 
targets  from  delay-doppler,  spotlight-mode  radar  data.  This  approach  is  based  on 
a  model  for  the  target’s  reflectivity  in  terms  of  wide-sense  stationary,  uncorrelated 
scatterers  having  complex-valued  Gaussian  statistics.  The  imaging  problem  is  to 
estimate  the  target’s  scattering  function  in  terms  of  radar-echo  data  acquired  with 
a  series  of  target  illuminations.  We  develop  [in  Appendix  6.3]  a  method  for  solving 
this  multidimensional  spectrum  estimation  problem  through  the  use  of  maximum- 
likelihood  estimation  implemented  via  the  expectation-maximization  algorithm.” 

2.2.4  Reprint  of:  P.  Moulin.  D.  L.  Snvder.  and  J.  A.  O’Sullivan:  "A  Sieve-Constrained 
Maximum-Likelihood  Estimator  for  the  Spectrum  of  a  Gaussian  Process".  Proc.  28th 
AHerton  Conference.  Urbana-Champaign.  IL.  Sept.  1989:  see  Appendix  6.4  for  a  reprint. 

"Maximum-likelihood  spectrum  estimation  is  an  ill-posed  problem.  In  this 
paper,  we  use  of  a  method  of  sieves  for  addressing  this  issue.  The  estimate  of  the 
spectrum  is  constrained  to  a  subset  of  some  Hilbert  space  of  functions  over  which 
a  complete  set  of  nonorthogonal  basis  functions  is  defined.  The  estimate  is  then 
represented  by  a  countable  set  of  coefficients  in  a  nonorthogonal  series  expansion. 

By  defining  an  appropriate  sieve  on  this  countable  set,  our  problem  reduces  to 
maximum-likelihood  estimation  of  the  parameters  in  the  sieve.  Three  main  attractive 
features  of  this  approach  are:  (1)  the  nonorthogonal  expansion  is  a  convenient 
framework  for  defining  the  sieve  and  including  a  priori  information;  (2)  mean-square 
consistency  of  the  estimates  can  be  expected;  and  (3)  we  have  derived  a  tractable 
alternating  maximization  algorithm  for  estimating  the  parameters.  The  setup  of  this 
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problem  is  general  and  can  be  applied  without  major  difficulties  to  the  estimation 
of  higher-dimensional  spectral  functions,  as  occurs,  for  example,  in  imaging  radar 
targets  from  delay-doppler  data." 

2.2.5  Abstract  of:  J.  A.  O’Sullivan.  P.  Moulin.  D.  L.  Snvder.  and  S.  P.  Jacobs.  "Computa¬ 
tional  Considerations  for  Maximum-Likelihood  Radar  Imagine."  Proc.  1990  CISS.  Prince¬ 
ton  University;  see  Appendix  6.5  for  a  reprint. 

"Recent  papers  have  outlined  a  new  approach  for  spectrum  estimation  and  radar 
imaging  based  on  expectation- maximization  algorithms  for  structured  covariance 
estimation.  Performance  of  this  approach  has  been  promising  for  the  problems 
studied.  Application  of  this  approach  to  real  data  sets  has  been  limited,  however, 
due  to  the  need  to  invert  a  matrix  whose  dimension  equals  the  size  of  the  data  set. 

For  radar  applications  where  an  image  is  to  be  formed,  data  sets  can  be  on  the  order 
of  2 14  for  128x  128  images.  This  makes  the  use  of  the  new  approach  difficult  in 
its  previously  described  form.  This  paper  proposes  both  approximation  methods  for 
inverting  typical  matrices  and  constraints  on  radar  transmitted  signals  which  make 
maximum  likelihood  image  estimation  viable.  These  constraints  may  be  satisfied 
for  real  signals  used  in  radar  imaging  systems.  Simulations  are  shown  to  demonstrate 
the  performance  of  the  algorithms.  Finally,  motivated  by  the  images  resulting  from 
the  simulations,  regularization  methods  are  discussed." 

3.  Other  Project  Activities 

3.1  Optical  Radar  Workshop 

At  the  request  of  Dr.  W.  J.  Miceli,  Office  of  Naval  Research,  Boston  MA,  we  organized 
and  hosted  a  one-day  workshop  on  laser  radar  imaging  on  April  12,  1989.  The  purpose  of  the 
meeting  was  to  discuss  various  tomographic  image  reconstruction  methods,  their  applicability 
to  laser  radar  imaging,  and  their  implementation  via  a  suitable  real-time  signal  processing 
system.  The  goal  was  to  provide  technical  interaction  among  researchers  interested  in  the  topic. 
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Most  participants  were  funded  through  O.N.R.  and/or  SDIO/T/IS.  The  approximately  twenty 
attendees  were  from  government  laboratories  and  organizations,  university  research  laboratories, 
and  industry. 

3.2  Invited  Presentations 

In  August  1989,  a  special  program  on  the  subject  of  signal  processing  was  held  at  the 
Institute  for  Mathematics  and  Its  Applications,  of  the  University  of  Minnesota  in  Minneapolis. 
We  were  invited  to  present  our  work  on  the  radar  imaging  problem  and  to  prepare  a  chapter 
for  a  book  to  be  published  by  the  I.  M.  A.  This  chapter,  coauthored  by  J.  O’Sullivan  and  D. 
Snyder  and  titled  "High  Resolution  Radar  Imaging  Using  Spectrum  Estimation  Methods,"  will 
appear  during  1990;  a  preprint  is  in  Appendix  6.6.  As  a  result  of  the  interest  in  radar  detection 
and  imaging  problems  that  developed  during  this  program,  a  second  program  on  the  topic  of 
radar  and  sonar  has  been  organized  and  will  take  place  in  June  1990  at  the  I.  M.  A. 

4.  Work  in  Progress 

4.1  Real-Data  Experiment 

An  effort  to  collect  real  data  with  which  to  test  and  compare  our  methods  for  radar  imaging 
has  been  initiated  in  collaboration  with  the  McDonnell-Douglas  Company  in  St.  Louis.  A 
sphere  having  a  diameter  of  one  meter  and  a  rough  surface  was  placed  on  a  rotating  pedestal 
in  a  compact  radar  test-range,  and  data  were  collected  at  several  values  of  signal-to-noise  ratio. 
This  object  was  selected  because  of  its  simplicity  and  the  fact  that  its  scattering  function  can 
be  predicted  analytically,  thereby  providing  a  test  object  having  known  characteristics  with 
which  to  compare  results.  These  data  have  only  recently  been  acquired  and  have  not  yet  been 
processed. 
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4.2  Graduate-Student  Theses 


4.2.1  Pierre  Moulin 

Pierre  Moulin  is  presently  writing  his  doctoral  dissertation  on  the  subject  of  estimation 
methods  for  forming  images  of  diffuse  radar  targets.  It  is  anticipated  that  this  thesis  will 
be  completed  in  May  1990. 

4.2.2  Kenneth  Krause 

Kenneth  Krause  is  presently  pursuing  research  for  his  doctoral  dissertation  on  the 
subject  of  forming  images  of  specular  radar  targets.  A  goal  is  to  develop  a  method  that 
accommodates  both  diffuse  and  specular  components  in  the  radar  echo.  It  is  anticipated 
that  this  thesis  will  be  completed  in  1991. 

4.2.3  Steven  Jacobs 

Steven  Jacobs  is  presently  writing  his  master’s  dissertation  on  the  subject  of  compu¬ 
tational  issues  associated  with  forming  radar  images  using  the  estimation  methods  we  have 
developed.  It  is  anticipated  that  this  thesis  will  be  completed  by  August  1990. 

5.  Personnel 

The  personnel  who  participated  in  this  research  project  during  the  reporting  period  are  the 
following. 

Steven  Jacobs 

•Graduate  Research  Assistant  in  the  Electronic  Systems  and  Signals  Research  Laboratory 

•  M.  Sc.  candidate  in  the  Department  of  Electrical  Engineering 

•  received  no  support  under  the  O.N.R.  Contract 

•task:  examine  computational  issues  in  radar  imaging  and  develop  parallel  implementation 
strategies 
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Kenneth  Krause 

-Part-Time  Graduate  Assistant  in  the  Electronic  Systems  and  Signals  Research  Laboratory 

-  Employed  by  the  McDonnell-Douglas  Astronautics  Co. 

•  Ph.  D.  Candidate  in  the  Department  of  Electrical  Engineering 

-  received  no  support  under  the  O.N.R.  Contract 

-  task:  include  specular  components  in  the  maximum-likelihood  method 
Pierre  A.  Moulin 

-Graduate  Research  Assistant  in  the  Electronic  Systems  and  Signals  Research  Laboratory 

-  Ph.  D.  Candidate  in  the  Department  of  Electrical  Engineering 

-  received  support  as  a  Graduate  Research  Assistant  under  the  O.N.R.  Contract 

-  task:  analyze  performance  of  the  maximum-likelihood  method,  include  sieve  and 
signal-to-noise  ratio  constraints  for  regularization,  examine  computational  issues  to  make 
the  method  practical 

Joseph  A.  O’Sullivan 

-  Faculty  Research  Associate  in  the  Electronic  Systems  and  Signals  Research  Laboratory 

-  Assistant  Professor  of  Electrical  Engineering 

-  received  support  as  a  Senior  Research  Associate  under  the  O.N.R.  Contract 

•  task:  participate  in  all  aspects  of  the  research  project 

Donald  Porter 

-  Undergraduate  Research  Assistant  in  the  Electronic  Systems  and  Signals  Research 
Laboratory 

-  B.S.E.E.  and  B.S.C.S.  degree  candidate  in  the  Electrical  Engineering  and  Computer 
Science  Departments 

•  received  no  support  under  the  O.N.R.  Contract 

•  task:  implement  conventional  radar  imaging  algorithms,  implement  a  computer  sim¬ 
ulation  of  radar  echo  data  based  on  a  given  scattering  function 


Donald  L.  Snyder 

•  Principal  Investigator 

•  Director  of  the  Electronic  Systems  and  Signals  Research  Laboratory 

-  Professor  of  Electrical  Engineering 

-  received  support  under  the  O.N.R.  Contract 

•  task:  participate  in  all  aspects  of  the  research  project 

Michael  Turmon 

•Graduate  Research  Assistant  in  the  Electronic  Systems  and  Signals  Research  Laboratory 

-  M.  Sc.  candidate  in  the  Department  of  Electrical  Engineering 

•  received  no  support  under  the  O.N.R.  Contract 

-  task:  study  implementation  of  spectrum-estimation  methods  of  use  in  radar  imaging 
on  a  massively  parallel  computer  architecture  (1024  element  A.M.T.  DAP  with  SUN/4 
host) 

J.  Trent  Wohlschlaeger 

•Graduate  Research  Assistant  in  the  Electronic  Systems  and  Signals  Research  Laboratory 

•  Ph.  D.  Candidate  in  the  Department  of  Electrical  Engineering 

-  received  no  support  under  the  O.N.R.  Contract  during  the  reporting  period 

•  task:  study  tomographic  methods  for  forming  radar  images 


-  to  - 


6.  Appendices 


6.1  Reprint  of:  D.  Snyder,  J.  O’Sullivan,  and  M.  Miller,  The  Use  of  Maximum  Likelihood 
Estimation  for  Forming  Images  of  Diffuse  Radar  Targets  from  Delay-Doppler  Data,  IEEE 
Trans,  on  Information  Theory ,  Vol.  35,  pp.  536-548,  1989. 


-  li  - 


T-IT/35/3//28182 


The  Use  of  Maximum  Likelihood 
Estimation  for  Forming  Images 
of  Diffuse  Radar  Targets  from 
Delay-Doppler  Data 


Donald  L.  Snyder 
Joseph  A.  O’Sullivan 
Michael  I.  Miller 


Reprinted  from 

1FFF.  TRANSACTIONS  ON  INFORMATION  THEORY 
VoL  35,  No.  3,  May  1989 


536 


[FEE  TRANSACTIONS  ON  INFORMATION  THEORY.  VOL.  35,  NO.  3,  MAY  1989 


The  Use  of  Maximum  Likelihood 
Estimation  for  Forming  Images 
of  Diffuse  Radar  Targets  from 
Delay- Doppler  Data 

DONALD  L.  SNYDER,  fellow,  ieee,  JOSEPH  A.  O’SULLIVAN,  member,  ifff., 

and  MICHAEL  I.  MILLER 


Abstract—  A  new  approach  to  high- resolution  radar  imaging  is  pre¬ 
sented.  The  starting  point  is  a  model  of  the  radar  echo  signal  based  on  the 
physics  governing  radar  reflections.  This  model  has  been  used  several 
times  in  the  past  for  describing  radar  targets  that  are  rough  compared  to 
the  wavelength  of  the  transmitted  radiation.  Without  specifying  precisely 
what  the  transmitted  signal  is,  a  general  estimation-based  procedure  is 
derived  for  obtaining  images.  After  discretizing  the  model,  the  radar 
imaging  problem  reduces  to  the  task  of  estimating  discretized  second-order 
statistics  of  the  reflectance  process  of  the  target.  Maximum  likelihood 
estimates  of  these  statistics  are  obtained  as  the  limit  point  of  an  expecta¬ 
tion-maximization  algorithm. 

I.  Introduction 

RADAR  SYSTEMS  can  be  used  to  produce  high-reso- 
lution  images  of  a  reflecting  target.  This  is  accom¬ 
plished  by  illuminating  the  target  with  a  series  of  pulses 
and  observing  the  return  echoes.  Each  patch  on  the  target 
introduces  a  certain  amount  of  propagation  delay  and 
Doppler  shift  to  a  pulse  it  reflects,  the  amount  depending 
on  the  range  and  velocity  of  the  patch  relative  to  the  radar 
transmitter  and  receiver.  The  beamwidth  of  the  radar 
antenna  relative  to  the  size  of  the  target  is  an  important 
consideration.  Images  can  be  produced  by  scanning  a 
narrowly  focused  beam  over  the  target  in  some  type  of 
raster  pattern  and  then  displaying  the  received  power  in 
delay  and  Doppler  or,  equivalently,  range  and  cross- range 
coordinates.  Images  can  also  be  formed  by  illuminating  the 
entire  target  in  spotlight  mode  with  a  wide,  relatively 
unfocused  beam.  The  received  signal  for  each  illumination 
is  then  a  complicated  superposition  of  the  echoes  received 
from  all  the  patches  that  make  up  the  extended  surface  of 
the  radar  target.  Our  concern  will  be  with  forming  images 
of  rotating  rough  targets  using  a  spotlight-mode  radar. 

We  will  denote  the  complex  envelope  of  the  signal 
transmitted  by  the  radar  by  (2£r)1/2rr(f),  where  £r  is  the 
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transmitted  energy,  and  sT(t)  is  normalized  to  unit  energy. 
The  particular  form  of  this  signal  will  not  need  to  be 
specified.  The  expressions  we  obtain  for  producing  an 
•  image  can  then  be  specialized  for  any  signal  of  choice, 
including  the  stepped-frequency  and  wide-band  chirp 
waveforms  used  in  practice,  as  discussed  by  Wehner  [1] 
and  Mensa  [2]  and  the  chirp-rate  modulated  waveforms 
discussed  by  Bemfeld  [3]  and  Snyder  et  al.  [4].1  When 
specializing  sT{t),  it  should  be  kept  in  mind  that  this 
represents  the  entire  sequence  of  transmitted  pulses  that 
illuminate  the  target. 

Walker  [5]  gives  a  clear  intuitive  description  of  the  radar 
imaging  problem.  He  considers  a  small  nonfluctuating 
reflector  rotating  counterclockwise  at  the  rate  of  fr  revolu¬ 
tions  per  second  on  a  circle  of  radius  r  centered  at  a 
distance  R0  from  the  radar  transmitter/receiver,  as  shown 
in  Fig.  1.  The  distance  to  the  reflector  at  time  t  is  given 


y,  cross  range 

i 


approximately  by 

R  ( t )  =  R0  +  x0cos  (2trfrt )  +  y0  sin  (2  vfrt ) , 

provided  R0  »  r  -  (x%  +  yj )1/2,  where  (x0,  y0)  is  the 
(x,  y)  position  of  the  reflector  at  time  t  =  0.  Then  the 
radar  echo  signal  sR(t)  received  following  an  illumination 
by  the  transmitted  signal  will  be  of  the  form 

J*(0  ~/2E^sT(t-r)b 

where  r»2 R(t)/c  is  the  two-way  propagation  delay  to 

xSote  added  in  proof  :  It  has  come  to  our  attention  since  submitting  this 
manuscript  that  chirp-rate  modulation  is  also  discussed  by  E  Feig  and  A. 
GrQnbaum  in,  “Tomographic  methods  in  range-Doppier  radar.”  Invent 
Probetms,  vol.  2.  pp.  185-195,  1986. 
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the  reflector,  with  c  being  the  propagation  velocity.  The 
quantity  b  is  a  complex-valued  scale  factor  which  models 
the  strength  of  the  received  echo.  This  scale  factor  will  be 
called  the  reflectivity.  It  includes  the  effects  of  inverse 
square-law  attenuation  experienced  by  the  propagating 
radiation  and,  importantly,  the  properties  of  the  reflector 
that  are  significant  in  the  electromagnetic  scattering  inter¬ 
action,  including  its  shape,  size,  and  surface  properties. 
More  generally,  the  reflectivity  can  vary  with  time  because 
the  aspect,  and  therefore,  the  scattering  interaction  with 
the  reflector  will  vary  as  it  rotates.  If,  as  discussed  by 
Walker  [5],  the  radar  data  sK(t)  are  examined  over  a  small 
interval  of  time,  then  the  delay  r  and  Doppler  shift  fD  can 
be  approximated  by 

r  “  2c-l(f?o  +  *o) 

and 


2  dR  2 


where  \  is  the  wavelength  at  the  carrier  frequency  of  the 
radar.  Thus  the  reflectivity  b ,  range  x0,  and  cross-range  y0, 
relative  to  the  coordinate  axis,  can  be  determined  from  the 
amplitude,  delay,  and  Doppler  information  contained  in 
the  radar  data.  Extracting  this  information  permits  the 
formation  of  an  image  of  the  reflector  by  displaying  |b|  or 
|b|2  at  the  appropriate  location  in  range  and  cross-range 
coordinates.  The  maximum  delay  and  Doppler  shift  are 
determined  by  the  distance  (x%  +  y£)1/2  of  the  reflector 
from  the  coordinate  center  about  which  it  rotates  and  the 
rotation  rate  /,;  more  generally,  the  extent  of  a  reflector  in 
delay  and  Doppler  is  determined  by  the  physical  extent  of 
the  reflector  and  the  rotation  rate. 

Now  consider  a  spatially  extended  target  that  is  rotating. 
A  patch  on  the  surface  with  a  two-way  delay  in  the  interval 
[t,  t  +  At)  reflects  a  signal  that  is  incident  on  the  patch  at 
time  t  with  a  reflectance  strength  b(t,  r)  At.  Consequently, 
the  complex  envelope  of  the  received  echo  signal  sR(t) 
following  the  illumination  of  the  target  by  sT(t)  is  given  by 
the  following  superposition  of  returns  from  reflecting 
patches  at  all  the  two-way  delays  r: 

2T’T)<iT-  M 


The  total  received  signal  r(t)  is  also  assumed  to  be  cor¬ 
rupted  by  an  independent  additive  noise 

',(0“**(0  +  w(0  (2) 


where  w(t)  is  a  complex- valued  white  Gaussian  process 
with  a  mean  of  zero  and  a  covariance  function  defined  by 

£[w(/)w*(/')]  =iV06(r-r')  (3) 

where  the  asterisk  denotes  complex  conjugation.  We  refer 
to  b{  t,  t)  as  the  reflectance  process.  This  is  a  complex- val¬ 
ued  random  process. 

There  are  two  images  that  may  be  displayed  as  the  result 
of  processing  r(t)  with  an  estimation  procedure.  One  is  an 
estimate  of  the  reflectance  process  b(t,r)  itself,  and  the 
other  is  an  estimate  of  the  covariance  or,  equivalently,  the 


spectral  density  of  this  process.  These  may  be  regarded  as 
conditional  first-  and  second-order  statistics  of  the  reflec¬ 
tivity  process,  respectively,  in  terms  of  the  radar  data  (2). 

If  b(t,  r)  is  deterministic,  define  c(f,  r)  to  be  its  Fourier 
transform  in  the  t  variable, 

C(/>T)=/  b(t,r)e~,2wf'  dt.  (4) 

The  function  c(f,r)  then  contains  information  about  the 
target  in  delay  t  and  Doppler  /  coordinates.  An  image  of 
the  target  is  obtained  by  placing  the  magnitude  or  squared 
magnitude  of  this  function  into  delay  and  Doppler  bins. 
We  refer  to  this  as  the  reflectance  image.  This  transform 
can  be  obtained  in  a  variety  of  ways,  depending  on  the 
signal  sT(t)  selected  to  illuminate  the  target.  For  the 
stepped- frequency  signals  used  in  practice,  the  usual  ap¬ 
proach  consists  of  two  operations  described  by  Wehner  [1]. 
The  first  is  to  place  the  data  into  delay  (or  range)  bins  by 
separately  Fourier  transforming  N  sample  values  of  the 
received  signal  acquired  for  each  transmitted  group  of  ,V 
stepped- frequency  pulses.  The  resulting  delay-binned  data 
are  placed  in  the  rows  of  an  N  x  N  matrix  where  each  row 
contains  the  transformed  data  from  one  pulse  group.  In 
the  second  operation,  the  columns  of  this  matrix  are 
Fourier  transformed  to  obtain  a  Doppler  (or  cross-range) 
profile  at  each  delay.  The  resulting  two-dimensional  array 
is  intended  to  be  a  discrete  version  of  c(f.r)  in  delay 
(range)  and  Doppler  (cross-range)  coordinates.  This  pro¬ 
cessing  based  on  two-dimensional  Fourier  transforms  is 
derived  using  a  strictly  deterministic  analysis  and  so  does 
not  account  for  statistical  properties  of  the  reflectance 
process  or  for  noise  that  may  be  present.  A  similar  process¬ 
ing  is  employed  for  the  linear  FM-chirp  signals  also  used 
in  practice  for  radar  imaging  [1],  [2], 

For  situations  in  which  the  target’s  surface  is  rough 
compared  to  the  wavelength  at  the  carrier  frequency.  b(t,r) 
may  be  taken  to  be  a  complex-valued  Gaussian  random 
process,  as  discussed  by  Shapiro  et  al.  [6]  for  radar  systems 
operating  at  laser  frequencies  and  Van  Trees  [7]  at  mi¬ 
crowave  frequencies.  If  there  are  no  glint  or  specular 
components  in  the  echo,  then  this  is  a  zero  mean  process 
with  covariance 

E[b{t.r)b*{t\r')\  =*  K{t-t',T)8{T-T').  (5) 

The  delta  function  in  this  expression  results  from  postulat¬ 
ing  that  each  reflecting  patch  introduces  an  uncorrelated 
contribution  to  the  echo  signal.  That  the  function 
K(t  -  t',  t)  depends  only  on  the  difference  of  t  and  t\  and 
not  on  t  and  t'  separately,  results  from  postulating  that 
the  reflectance  process  is  wide-sense  stationary  for  each 
delay.  A  reflectance  process  with  these  properties  is  said 
by  Van  Trees  [7]  to  possess  wide-sense  stationary  uncorre- 
lated  scatterers  (WSSUS).  Assuming  that  the  reflectance 
process  has  these  properties,  the  delay-Doppler  data  asso¬ 
ciated  with  the  radar  target  may  be  obtained  from  the 
Fourier  transform  of  K(t,  t)  in  the  t  variable, 

S(/.t)«  (**  K(t,r)exp(  -  jlirft )  dt.  (6) 

J  -  ao 
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The  function  S(/,  t)  is  called  the  scattering  function  of 
the  target  and,  as  a  function  of  /,  is  the  power-density 
spectrum  of  the  reflectance  process  at  each  delay  r. 
S(f.  t)  A/At  is  the  mean-square  strength  (or  power)  of  the 
reflectance  of  all  patches  on  the  target  having  a  Doppler 
shift  in  the  interval  [/,  /  +  A/)  and  a  delay  in  the  interval 
[t,  t  +  At).  The  scattering  function  may  be  viewed  in  delay 
and  Doppler  coordinates  as  an  image  of  the  target.  We  call 
this  the  scattering  function  image. 

Our  approach  to  forming  radar  images  will  be  to  use 
maximum  likelihood  methods  with  the  data  model  in  (2)  to 
estimate  the  scattering  function.  We  will  also  obtain  an 
estimate  of  the  reflectance  process.  Model-based  ap¬ 
proaches  that  use  statistical  estimation  theory  techniques 
to  derive  image-formation  algorithms  appear  less  fre¬ 
quently  in  the  literature  about  radar  imaging  than  the 
deterministic  approaches  outlined  above.  One  example  is 
that  of  Frost  et  al.  [8],  who  uses  a  multiplicative  model  and 
Wiener  filtering  techniques.  The  approach  we  will  describe 
differs  in  that  the  model  (2)  we  adopt  of  the  echo  signal  is 
more  complicated  than  a  simple  multiplicative  one  and 
depends  explicitly  on  the  transmitted  waveform  through  a 
spatial  integration  over  the  reflecting  target.  We  also  do 
not  restrict  the  processing  to  be  linear;  in  particular,  we 
show  that  the  processing  of  the  received  data  for  produc¬ 
ing  the  maximum  likelihood  estimate  of  the  scattering 
function  and  a  corresponding  estimate  of  the  reflectance 
process  is  not  linear.  An  approach  for  estimating  scattering 
functions  of  spread  channels  is  given  by  Gaarder  [9],  who 
cites  earlier  work  on  the  subject  by  Green  [10],  Kailath 
[11],  Gallager  [12],  Hagfors  [13],  [14],  Price  [15],  Levin  [16], 
Abraham  [17],  Sifford  [18],  and  Reiffen  [19].  Gaarder 
assumes  a  specific  processing  architecture  in  the  form  of  a 
cascade  of  a  linear  filter  square-law  envelope  detector  and 
another  linear  filter  and  claims  that  this  processing  is 
either  more  general  than  or  equivalent  to  those  of  most 
previous  authors.  Our  approach  differs  in  that  no  particu¬ 
lar  processing  is  assumed  in  advance;  rather,  we  derive  the 
processing  to  produce  the  estimates,  starting  from  a  model 
for  the  received  data  and  applying  recent  results  in  maxi¬ 
mum  likelihood  estimation.  The  processing  which  results  is 
quite  distinct  from  that  discussed  by  Gaarder. 

For  our  new  approach  to  radar  imaging,  we  adopt  the 
WSSUS  model  of  a  diffuse  radar  target  described  by 
Shapiro  et  al.  [6]  and  Van  Trees  [7].  We  treat  both  the 
reflectance  process  and  its  second-order  statistic,  the  scat¬ 
tering  function,  as  unknown  quantities.  The  iterative  ap¬ 
proach  we  develop  for  forming  images  yields  the  maxi¬ 
mum-likelihood  estimate  of  the  scattering  function  and, 
simultaneously,  the  conditional-mean  estimate  of  the  re¬ 
flectance  process  based  on  statistics  which  are  consistent 
with  the  estimated  scattering  function.  Thus  both  of  the 
quantities  treated  separately  in  other  imaging  schemes  are 
produced  simultaneously  with  our  new  approach,  which  is 
a  unique  and  important  aspect  of  our  approach. 

We  will  develop  a  necessary  condition,  called  the  trace 
condition,  which  the  maximum  likelihood  estimate  of  the 
target’s  scattering  function  must  satisfy.  This  equation 


appears  to  be  very  hard  to  solve  analytically.  As  a  conse¬ 
quence,  we  reformulate  the  imaging  problem  using  the 
concept  of  incomplete-complete  data  spaces  and  then  use 
the  expectation-maximization  algorithm  of  Dempster  et  al. 
[20]  to  derive  an  iterative  algorithm  for  producing  the 
maximum  likelihood  estimate  of  the  scattering  function. 
The  technique  we  use  to  accomplish  this  parallels  that 
described  by  Miller  and  Snyder  [21]  for  power-spectrum 
estimation  and  extends  their  work  to  include  indirect  mea¬ 
surements  of  the  process  whose  spectrum  is  sought;  the 
process  is  now  measured  following  the  linear  transforma¬ 
tion  and  additive  noise  seen  in  (2).  As  shown  by  Turmon 
and  Miller  [22],  this  is  a  high-resolution  approach  to  spec¬ 
trum  estimation  which  results  in  estimated  spectra  with 
smaller  bias  and  mean-square  error  than  other  recently 
developed  approaches  discussed  in  the  literature.  We  ex¬ 
pect  that  similar  improvements  will  be  seen  in  radar  im¬ 
ages  of  scintillating,  diffuse  targets  when  this  new  tech¬ 
nique  is  used. 

II.  Maximum  Likelihood  Imaging  for  the 
Incomplete  Data  Model 

For  reasons  that  will  become  evident  in  the  next  section, 
we  term  the  data  r(t)  in  (2)  the  incomplete  data  for  the 
radar-imaging  problem.  The  model  given  in  the  Introduc¬ 
tion  for  these  data  consists  of  the  sum  of  the  radar  echo 
signal  sR(t )  of  (1)  and  an  independent  white  noise  process 
w(i).  We  may,  therefore,  state  the  problem  of  imaging  a 
diffuse  radar  target  as  that  of  estimating  the  scattering 
function  S(f,r)  or  equivalently,  the  covariance  function 
K{t,  t)  given  radar-return  data  { r(t ),  T,£tzTf}  on  an 
observation  interval  ( TnTf ).  In  this  section,  we  first  dis¬ 
cretize  the  model  for  the  incomplete  data  and  then  derive 
and  discuss  a  necessary  condition,  called  the  trace  condi¬ 
tion,  which  the  maximum  likelihood  estimate  of  the  dis¬ 
cretized  version  of  K(t,  r )  must  satisfy. 

Discrete  Model 

In  anticipation  of  using  discrete-time  processing  of  radar 
data  to  produce  images,  we  now  state  the  discrete  version 
of  our  model  as  follows.  We  are  given  N  samples  of  the 
complex-valued  radar  data  corresponding  to  (2), 

r(n)  =  sR(n)+ w(n),  n  =  0,1,- •  • ,  N  —  1  (7) 

where  w(n)  is  a  white  Gaussian  sequence  with  zero  mean 
and  covariance 

£[w(n)w*(n')]  (8) 

where  S„  is  the  Kronecker  delta  function,  and  the  signal 
samples  corresponding  to  (1)  are  given  by 

_  +  30 

sR(n)-J2ET  £  sT(n,i)b(n,i), 

I  —  “  » 

n-0,l.-”,JV-l.  (9) 

In  this  expression,  we  define  sT(n.  i)  and  b(n.  i)  in  terms 
of  the  transmitted  signal  and  the  reflectance  process,  re- 
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where  the  iV-dimensional  vectors  sR  and  w  are  given  by 

s*( 0)  w(0)  ' 

w(0  ,  .. 

sR=  .  w=  (15) 


Sr(Q) 

AO) 

**(1) 

w  = 

w(l) 

[AN- 1), 

spectively,  according  to 

sT(n ,i)  =  sT(n  Ar  -  i  At) 


b(n,  i)  =  6  n  Af  -  —i  At,  i  At  At  (11) 

\  2  j  Also,  define  5*  as  the  NIR  x  N  rectangular  matnx  ex¬ 

pressed  in  column-block  form  in  terms  of  IR  separate 
N  X  N  matrices  according  to 

where  A i  and  At  are  the  sampling  intervals  adopted  m  the 
discretization  in  time  and  delay,  respectively.  We  assume 

that  the  target  has  a  finite  extent  in  range;  thus  b(rt,i)  and  Si 

therefore  terms  forming  the  sum  in  (9)  are  zero  for  i  ^  =  :  ' 

outside  the  /„  (here,  the  subscript  R  denotes  range)  values  ^  ' 

m,  m  +  1,  •••,«  +  IR -1  starting  from  the  minimum  two-  r*~1' 

way  delay  corresponding  to  m.  This  discrete  reflectance  is  where  Sj  is  an  N  x  N  diagonal  matrix  containing  sample 

a  Gaussian  sequence  with  zero  mean  and  covariance  given  values  of  the  complex  envelope  of  the  transmitted  signal 

— -  sT(t), 


sr(0,  m  +  j) 

0  0  • 

0 

Q 

jt(1,  m  +  j)  0  • 

0 

0 

0  •  • 

0 

0 

0  •  • 

•  sT(N  —  l,  m  +  j) 

£[6(n,i)6*(n',i')]  -  K(n  -  n',  i)  8, ,  r.  (12) 

The  discrete  scattering  function  S(v,  i)  is  the  Fourier 
transform  of  K(n,i)  in  the  n  variable, 

+  00 

S(v,i)  =  Y,  K(n,i)exp(- jlnvn).  (13) 


The  imaging  problem  for  the  discrete  model  is  to  estimate 
S(u, /)  or,  equivalendy,  the  covariance  function  K(n,  i), 
for  all  frequencies  o  spanning  the  target  in  Doppler,  and 


Further,  define  the  reflectance  vector  b  of  dimension  ;V/R 
in  the  column-block  form  of  IR  vectors  according  to 

m 

6(1) 

6-  V  (18) 

where  each  6(i)  is  a  vector  of  dimension  N, 

6(0,  m  +  i) 

,  6(1,  m  +  i) 

MO*  ;  •  (19) 

6(JV-l,m-H) 


for  all  delays  i  spanning  the  target  in  the  delay,  given  the  Using  (9)  and  these  definitions,  we  can  now  express  the 
radar  data  {/•(«),  n  »  0,1,- •  *,  N  -1}.  Af-dimensional  signal  vector  sR  of  (14)  and  (15)  as 


Matrix  Model 


sR  =  j2E^S+b 


These  discrete  equations  may  conveniently  be  written  in  w^ere  a  superscript  plus  sign  denotes  the  Hermitian-trans- 
matrix  form  as  follows.  Define  r  to  be  the  received-signal  Posc  °Pcration.  In  terms  of  these  defined  matrices,  the 
vector  of  dimension  N,  received  vector  has  zero  mean  and  covariance 

r(0)  Kr- E{rr+)-E{sRsZ)+  E{ww+) 

r(  1)  ^  —  2ErS* E(bb*  )S  +  N0I.  (21) 

r  ™  "  s*  w  (  ^  Then,  since 

r(N  —  1)  E(b(i)b*  (j))  -  K(i)  5,  y  (22) 

_ _  where  K{i)  is  the  Hermitian-symmetric  Toeplitz  matrix 

I  K(0,m  +  i)  K*(l,m  +  i)  •••  **(  N- 1,  m  +  i)  \ 


1  s.  +  w 


AN- 1), 


K(0,m  +  i) 
K(l,m  +  i) 

K(N  —  \,m  +  i) 


K*(\,m  +  i) 
K( 0,  m  +  i) 


K(0,m  +  i) 
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it  follows  from  (21)  that  the  covariance  Kr  of  r  is  given  by 
Kr  —  2EtS*KS  +  N0l  (24) 

where  K  is  the  block-diagonal  NIR  x  A'/^-dimensional  ma¬ 
trix  defined  by 


'K(O) 

0 

0  • 

0 

0 

K(  1) 

0  ■ 

0 

•  (25) 

I  0 

0 

0  • 

\  A  fA 

•  *(/*-!), 

The  /  th  diagonal  block  K(i)  of  K  is  the  covariance  matrix 
of  the  reflectance  process  at  the  ith  delay  bin. 

The  Estimation  Problem 

We  will  use  the  following  definition. 

Definition:  Let  K  denote  the  set  of  all  NIR  x  NIR 
block-diagonal  matrices  (25)  with  each  block  K(i)  an 
N  x  N  Hermitian-symmetric  Topelitz  matrix  (23).  Let  12  C 
K  be  a  specified  convex  subset  of  K.  Any  matrix  K  €  12  is 
termed  admissible.  A  variational  matrix  SK  €  K  is  called 
an  admissible  variation  of  K  for  a  fixed  K  e  12  if  there 
exists  an  a  >  0  such  that  K  +  SK  e  12  for  all  /J  satisfying 
0  <  /?  <  a. 

The  problem  is  to  form  an  admissible  estimate  of  the 
covariance  matrix  K  of  (25)  given  the  data  vector  r  of  (14). 
The  radar  image  then  viewed  is  the  discrete  scattering 
function,  an  estimate  of  which  may  be  obtained  from  the 
estimate  of  AT  by  use  of  (13). 

In  the  foregoing  definition  the  constraint  that  K  be  in  12 
is  used  to  obtain  a  “reasonable”  setup  of  the  problem. 
Here  we  assume  that  the  scattering  function  S(f,  r)  in  (6) 
is  only  nonzero  for  frequencies  /  satisfying  |/|  ^  for 

some  finite  upper  frequency  f _ and  for  all  delays  i.  This 

is  equivalent  to  the  assumption  of  a  target  of  finite  cross 
section  and  rotation  rate.  The  discrete-time  scattering 
function  S(v,i)  of  (13)  is  then  a  periodic  function  of  u 
consisting  of  a  sum  of  shifted  replicas  of  S(f,  r)  scaled  in 
amplitude  by  1/A/  and  in  frequency  by  A/,  where  A/  is  the 
time  between  samples  of  r(t).  The  replicas  of  S{f,  r)  are 
centered  at  every  integer  on  the  o  scale.  To  guarantee  that 
there  is  no  aliasing,  assume  that  the  sample  rate  1/A/  of 
r(t)  satisfies  the  Nyquist  condition  1/A/  >  2/^.  Then 
S(u,  /')  will  be  nonzero  between  - 1/2  and  + 1/2  only  for 
u  satisfying  |o|  ^  At.  The  output  of  our  algo¬ 

rithm  is  S(o,  i )  discretized  in  frequency  o.  For  a  resolution 
having  at  least  ^CR  (here,  the  subscript  CR  denotes  cross 
range)  samples  in  the  frequency  range  -  £  o  £  v^,  a 

total  of 

Am 

2  A  //_ 

samples  of  u  between  - 1/2  and  + 1/2  are  required. 

The  model  of  the  incomplete  data  r  of  (14)  is  that  r  is 
normally  distributed  with  zero  mean  and  covariance  speci¬ 
fied  in  (24).  Given  the  incomplete  data,  we  wish  to  esti¬ 
mate  the  covariance  K  of  the  reflectance  process,  as  de¬ 
fined  in  (25).  To  do  this,  we  adopt  the  maximum  likelihood 


method  of  statistics,  which  selects  K  to  maximize  the 
incomplete  data  log  likelihood 

Lid(  K )  =  -In(det(2£7-S+/fS  +  N0I)) 

-r+(2ETS+KS  +  N0l)~lr  (26) 

where  the  maximization  is  subject  to  the  constraint  that  K 
be  an  admissible  matrix. 

Lemma  1:  A  necessary  condition  for  an  admissible  K 
in  the  interior  of  8  to  be  a  local  maximum  of  Lid(K)  over 
all  K  s  8  is 

lt((2ETS+KS  +  N0I)~\rr*  -2ETS+ KS  -  N0I) 

■(2ErS+KS  +  NoI)~lS+SKS)=0  (27) 

for  all  admissible  variations  SK. 

The  proof  of  Lemma  1  in  the  Appendix  is  based  on  the 
fact  that  the  necessary  condition  for  an  admissible  K  to 
maximize  Lid(K)  is  that,  for  all  admissible  variations  SK, 


lim 

a  10* 


Lii(K  +  a8K)-Lid(K) 


(28) 


We  call  (27)  the  trace  condition.  Burg  et  al.  [23]  have 
studied  an  equivalent  problem  of  Toeplitz-constrained 
covariance  estimation  and  have  derived  the  trace  condition 
using  a  different  approach. 

If  C2  *  K ,  there  are  NIR  unknowns  in  K.  Since  SK  e  K, 
there  are  NIR  parameters  in  SK  that  can  be  varied  for  this 
case.  These  variations  in  the  trace  condition  generate  NIR 
equations  in  the  unknown  elements  of  K.  Thus,  in  princi¬ 
ple,  the  trace  condition  produces  enough  equations  to 
determine  the  unconstrained  maximum  likelihood  estimate 
K.  However,  the  equations  are  complicated  due  to  the 
inverse  matrices  appearing  in  (27);  thus  it  does  not  appear 
feasible  to  determine  K  directly  from  the  trace  condition. 
This  motivates  the  development  of  the  iterative  approach 
presented  in  the  next  section.  A  sequence  of  estimates  that 
increase  the  likelihood  at  each  iteration  stage  is  identified, 
and  we  demonstrate  that  stable  points  of  the  iteration 
satisfy  the  trace  condition. 

The  trace  condition  is  only  a  necessary  condition  that 
the  estimate  K  must  satisfy.  For  it  to  be  sufficient  as  well, 
the  second  derivative  must  be  negative  along  all  admissible 
variational  directions  SK. 

Lemma  2:  Sufficient  conditions  for  an  admissible  ma¬ 
trix  AT  in  the  interior  of  Q  to  be  a  local  maximum  of 
Lid(Af)  are  that,  first,  the  trace  condition  (27)  is  satisfied 
for  all  admissible  variations  SK  and.  second,  that 

tr  ( K;  'S  *  sksk;  1  ( 2  EtS*KS 

+  NqI  -  2rr*  )  Kf  lS*  SKS )  <  0  (29) 

for  all  admissible  variations  SK. 

The  proof  of  Lemma  2  is  given  in  the  Appendix.  Equa¬ 
tion  (29)  is  just  the  second  derivative  of  Ltd(K)  in  the 
direction  SK. 
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III.  Maximum  Likelihood  Imaging  for  the 
Incomplete/Complete  Data  Model 

The  fact  that  the  trace  condition  (27)  cannot  be  solved 
directly  for  the  maximum  likelihood  estimate  of  K  moti¬ 
vates  the  indirect  approach  we  now  take  of  embedding  the 
imaging  problem  in  a  larger,  seemingly  more  difficult 
problem.  The  result  will  be  an  iterative  algorithm  which, 
when  implemented,  produces  a  sequence  of  admissible 
matrices  /C<0\  K(1\- •  K(k),  •  •  ■  with  the  property  that 

the  corresponding  sequence  of  incomplete  data  log  likeli¬ 
hoods  Lid[tf<°>],  is  nonde¬ 

creasing  at  each  stage. 

Fuhrmann  and  Miller  [24]  have  recently  shown  that 
maximum  likelihood  estimates  of  Toeplitz-constrained  co- 
variances  which  are  positive  definite  do  not  always  exist 
when  given  only  one  data  vector  r.  A  necessary  and 
sufficient  condition  for  the  likelihood  function  to  be  un¬ 
bounded,  and  therefore  for  no  maximum  likelihood  esti¬ 
mate  to  exist,  is  that  there  be  a  singular  Toeplitz  matrix 
with  the  data  in  its  range  space.  For  our  imaging  problem, 
this  condition  is  that  an  admissible  K  exists  with 

2 ETS+kS+N0I  ' 

singular  so  that 

r  =  (2ETS+KS  +  N0I)a  (30) 

for  some  complex-valued  vector  a.  In  fact,  with  only  a 
Toeplitz  constraint  on  K,  a  sufficient  condition  that  a 
singular  estimate  for  K  be  obtained  is  that  N0  «  0  and  that 
a  singular  K  exist  with  r  in  the  range  space  of  2 ETS* KS. 
The  argument  for  this  mirrors  that  of  Fuhrmann  and 
Miller  in  [24,  theorem  1]  but  is  applied  to  the  complete 
data  log  likelihood  given  in  (A7)  of  the  Appendix. 
Fuhrmann  and  Miller  also  showed  that,  even  if  the  true 
covariance  had  eigenvalues  bounded  from  above  and  be¬ 
low,  the  probability  that  a  singular  Toeplitz  matrix  exists 
with  the  data  in  its  range  can  be  very  close  to  one.  By 
restricting  the  search  to  Toeplitz  matrices  with  circulant 
extensions,  they  were  able  to  show  that  the  probability  a 
singular  circulant  Toeplitz  matrix  has  the  data  in  its  range 
space  is  zero.  Thus,  for  maximum  likelihood  estimates  to 
be  nonsingular  with  probability  one  for  all  nonnegative 
values  of  N0 ,  we  may  restrict  the  class  of  admissible 
Toeplitz  matrices  to  be  those  with  circulant  extensions  of 
period  P,  where  P  is  equal  to  or  greater  than  the  number 
IV  of  data  samples  available  P>  N.  This  is  not  a  severe 
restriction  because  the  set  of  all  Toeplitz  matrices  is  ap¬ 
proached  by  the  subset  having  circulant  extensions  of 
period  P  as  P  tends  to  infinity.  What  we  envision  in 
adopting  this  constraint  is  that  for  each  delay  i,  the  N 
sample  values  of  the  reflectance  b(n,i),  n  «  0, 1,-  •  ■,  N  - 1, 
are  from  a  stationary  process  that  is  periodic  with  period 
P,  where  P  could  be  some  large  but  finite  value;  a  lower 
bound  on  P  in  terms  of  a  desired  cross-range  resolution  is 
discussed  above.  These  N  sample  values  of  the  reflectance 
enter  the  incomplete  data  r  according  to  (14)  and  (20).  By 
using  the  expectation-maximization  algorithm  of  Dempster 
et  at.  [20],  we  shall  develop  a  sequence  of  admissible 


matrices  that  have  the  maximum  likelihood  estimate  of  K 
subject  to  this  circulant  extension  as  a  stable  point.  The 
approach  parallels  that  of  Miller  and  Snyder  [21]  for 
estimating  the  power  spectrum  of  a  time  series  from  a 
single  set  of  data.  An  important  benefit  of  introducing  the 
periodic  extension  and  using  the  expectation-maximization 
algorithm  is  that  estimates  of  both  the  scattering  function 
and  the  reflectance  processes  are  obtained  simultaneously 
and  can  be  readily  viewed  as  target  images  in  range  and 
cross-range  coordinates;  thus  the  procedure  proposed  may 
be  considered  natural  for  the  imaging  problem  because 
both  types  of  images  considered  separately  in  the  past  are 
obtained  directly.  As  a  final  comment  regarding  our  use  of 
a  circulant  extension  for  K ,  we  note  that  in  estimating  a 
discretized  version  of  the  target’s  scattering  function,  the 
class  of  admissible  K  is  restricted  automatically  to  consist 
of  those  with  circulant  extensions.  For  completeness,  we 
also  include  in  the  Appendix  the  equations  obtained  using 
the  expectation-maximization  algorithm  for  estimating 
general  Toeplitz  matrices  when  the  assumption  of  a  circu¬ 
lant  extension  is  not  made. 

We  shall  introduce  a  modification  of  our  notation  to 
indicate  that  the  N  samples  of  the  reflectance  process  are 
from  a  stationary  periodic  process  of  period  P.  To  this  end 
let  i,v(0  denote  the  jV-dimensional  vector  b(i)  of  (19).  We 
now  think  of  f>v(i)  as  an  jV-dimensional  subvector  of  the 
/’-dimensional  vector  bP(i)  formed  from  samples  of  the 
reflectance  process  over  a  full  period. 

b(0,m  +  i) 

b  ( 1 ,  m  +  / ) 


MO 


b(  N  -  1,  m  +  i ) 


(31) 


[b(P -\,m  +  i)J 

If  Is  is  the  N  X  N  identity  matrix,  and  JR  the  P  X  ,V 
matrix  defined  by 


then 


(32) 


MO  =  A  bP(i). 

Also,  let  bv  denote  the  W-dimensional  vector  b  of  (18), 
and  bP  the  /’/-dimensional  vector  with  /  th  block  element 
bp(i).  Then, 


b„-M;bp 


where  MR  is  the  PIR  x  SIR  block-diagonal  matrix 


M 


K  = 


(JR  0  0 

0  JR  0 


0  ' 
0 


(33) 


\  0  0  0-  JR  j 

Furthermore,  let  K„(i)  denote  the  /V  x  N  Toeplitz  covari¬ 
ance  matrix  K{i)  of  bv(i)  defined  in  (23),  and  KR(i)  the 
P  X  P  circulant  covariance  matrix  of  bp(i).  Then,  the 
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Toeplitz  matrix  ATv(i)  is  the  upper  left  block  of  the 
circulant  matrix  KP(i),  as  given  by 

*„(«)  = MOV 

Lastly,  let  Kp  denote  the  PIR  x  PIR  block-diagonal  matrix 
in  the  form  of  (25)  with  the  /th  diagonal  block  being 
KP(i).  Then,  if  KN  denotes  the  NIRXNIR  matrix  K  of 
(25), 

k„=m;kpmr. 

Let  W  denote  the  P  x  P  discrete  Fourier-transform  ma¬ 
trix  scaled  so  that  the  columns  are  orthonormal, 


f*V° 

• 

o 

"p 

k  2k 

Wp  Wp 

tv/- 1)4 

w° 

\Wp 

h-;-1  • 

• 

M(P-lXP-i) 

Wp 

where  wp  =  exp(  -  jlir/P).  Also,  let  Wp  be  the  PI R  X  PIR 
block-diagonal  matrix 

W  0  0  •  •  0 

wP  =  0  w  0  ;  ;  0  .  (35) 

t  0  0  0  •  •  w, 

Then,  bP  can  be  represented  in  rotated  coordinates  accord¬ 
ing  to 

(  a(0)  \ 


where  a(i)  =  IVbp(i).  The  assumption  that  bP(i)  origi¬ 
nates  from  a  periodic  process  implies  that  the  P/R-dimen- 
sional  vector  aP  is  normally  distributed  with  zero  mean 
and  diagonalized  covariance 

AP  =  E{aPa+P)  =  WPKPW; .  (37) 

We  will  denote  the  (p  +  iP) th  diagonal  element  of  APby 
a2(i)\  this  is  the  p th  diagonal  element  of  the  PxP 
diagonal  matrix  E[a(i) a  *(/)). 

Let  S(u,i)  be  discretized  in  frequency  with  P  samples 
taken  for  0  ^  v  <  1.  These  samples  may  be  obtained  from 
(13)  as 

$(£,«)-  Oexp  (  -  j—~p~ )  (38) 

for  p  -  0, 1,  •••,/*- 1.  The  pth  such  sample  is  just  the  /Hh 
entry  in  the  vector 

/ PWKP(i)e  (39) 

where  e  is  the  /’-dimensional  unit  vector 


Substituting  (37)  into  (39),  we  get 

{P  WKf(i)e  =  {PA{i)We,  (40) 

but  Pl/2We  is  a  /’-dimensional  vector  of  ones,  and  there¬ 
fore 

s(£.‘)-«,2(0.  (41) 

which,  according  to  the  above  definition,  is  the  ( p  +  iP) th 
diagonal  element  of  the  diagonal  matrix  AP.  The  entries  of 
the  diagonal  matrix  Ap  in  (37)  are  then  samples  of  the 
scattering  function. 

The  constraint  from  Section  II  that  the  scattering  func¬ 
tion  S(o,.i )  be  nonzero  only  for  |v|  ^  =  ym«>  f°r 

values  of  v  between  -1/2  and  +1/2,  may  be  incorpo¬ 
rated  at  this  point  in  the  development.  Since  A,,  is  a 
diagonal  matrix  of  samples  of  the  scattering  function,  we 
restrict  AP  to  have  nonzero  values  only  in  its  top  left  and 
bottom  right  comers.  More  precisely,  let  /CR  be  the  small¬ 
est  odd  integer  satisfying 

/cr  >  2vtMX/>- 

/CR  is  the  number  of  cross-range  resolution  ceils  implied 
by  P  and  Then,  let  /cr  be  the  /cr  *  /*  matrix 

/  f/l  °  °\ 
lo  0  /J 

where  /j  is  an  [(/CR  +  l)/2)x[(/CR  +  l)/2)  identity  matrix 
and  I2  is  an  [(/CR  -  l)/2]x[(/CR  - 1)/2]  identity  matrix. 
Let  AfCR  be  the  ICKIR  x  PIR  matrix 

yCR  0  0  •  •  0 

Mck  »  0  /cr  0  1  •  0 

o  o  o  yCR 

Define  1P  to  be  the  diagonal  matrix 

^■p  ~  ^cr ApMqk- 

The  diagonal  elements  of  2  p  are  the  potentially  nonzero 
diagonal  elements  of  Ap.  Recognizing  that  some  dements 
of  the  diagonal  matrix  AP  are  zero  and  using  the  defini¬ 
tion  of  we  conclude  also  that 

Ap  =  M£K1PMCK. 

The  set  £2  referred  to  in  the  definition  in  Section  II  can 
now  be  specified.  We  restrict  consideration  to  those  covari¬ 
ance  matrices  generated  by  2P  from  above,  so 

Q  -  { K  e  K:  K  -  M;  W;  M^2pMckWpMr } . 

For  use  with  the  expectation-maximization  algorithm, 
we  identify  the  complete  data  as  ( cp ,  tv),  where  tv  is  the 
//-dimensional  noise  vector  defined  in  (15)  and  cP  is 
defined  by 

cp  ”  3/cr ap- 

Since  elements  of  ap  corresponding  to  the  zero  diagonal 
elements  of  A  p  are  almost  surely  zero,  we  see  also  that 

ap  ”  cp. 


13  - 
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Using  this  fact,  we  note  from  (14)  and  (20)  that  the 
incomplete  data  r  can  be  obtained  from  the  complete  data 
according  to  the  mapping 

r=*Y*cP+w  (42) 

where  we  define  the  IKICR  x  N  matrix  (which  appears 
throughout  the  development  that  follows): 

r  =  {Te^.mckwpmrs. 

The  log  likelihood  function  L^CZp)  of  the  complete  data 
as  a  function  of  2,,,  the  diagonal  covariance  matrix  of  cP, 
is  given  by 

Lca(2P)  =  —  ln(det  (2/,))  —  cp^.plcp 
/,-i  'cR-t 

”-2  £  £  ln(<J?(i)) 

i  —  0  p  —  0 


Since  Lcd  in  (43)  is  a  function  of  the  nonrandom  matrix 
2,,,  the  result  of  this  conditional  expectation  is  a  function 
of  2P.  It  is  also  a  function  of  2)>A,  because  the  expectation 
is  performed  assuming  that  cp  is  normally  distributed  with 
covariance  1{P\  We  have  indicated  this  dependency  of  the 
conditional  expectation  on  both  and  2),*’  in  the  defini¬ 
tion  of  the  function  Q  in  (44).  From  (43).  we  have  that 

0[2,|2J,*>]--2  £  £  !«(*,(/)) 

i  —  0  p  —  0 

i-0  p -  0 

(45) 

The  M-step  yields  the  estimate  I.{P+l)  at  stage  k  + 1  as  the 
choice  of  "2P  that  maximizes  this  conditional  expectation. 


-  £  £  M0lV2(0  (43) 

i-0  p—0 

where  all  terms  that  are  not  a  function  of  2P  have  been 
suppressed  and  cp(i)  is  the  p\h  element  of  the  /CR-dimen- 
sional  vector  cP(i)  =  JCKa(i). 

The  expectation-maximization  algorithm  for  estimating 
the  covariance  of  the  reflectance  process  KP  from  the 
incomplete  data  r  is  an  alternating  maximization  proce¬ 
dure  in  which  a  sequence  of  estimates  2)?\  2J}\*  •  •, 
2(/\  ■  -  •  of  1P  is  obtained  First,  where  the  expectation- 
maximization  procedure  specifies  how  to  obtain  2(/>*  +  l) 
from  2{pk)  for  k  =  0, 1,  •  •  • .  If  2),*'  denotes  the  estimate  of 
1P  at  stage  k,  then  there  is  a  corresponding  term 

K(pk)  -  w;M^ipMCPwP 

in  a  sequence  of  estimates  of  Kp.  Likewise,  to  the  kth  term 
KPk)  of  the  sequence  of  estimates  of  KP,  there  is  a  corre¬ 
sponding  term 

=  M+RKpMp _ 

Kj,k+"{ 0)  0 

o  x:<*+1,(i) 
o  o 


2(/  +  1>  =  argmax[2(2/,|2(p))].  (46) 

subject  to  the  constraint  that  the  maximizer  be  a  diagonal 
covariance  matrix.  From  (45),  this  maximization  yields  the 
diagonal  matrix  2(/>*  +  1)  with  ( p  +  //CR)th  diagonal  element 

(°/(0f +,>  =  E[|c,(,)|V.2^].  (47) 

Thus,  we  may  write  2{P* 11  as 

2</*1'=E[c,c;ir.2<;»]  (48) 

where  the  d  over  the  equal  sign  means  that  the  diagonal 
terms  in  the  matrix  on  the  left  side  equal  the  diagonal 
terms  in  the  matrix  on  the  right  side  and  that  all  the  off 
diagonal  elements  on  the  left  side  are  zero. 

Expression  (48)  appears  to  be  complicated  because  of 
the  several  matrices  we  have  defined,  but  it  produces  a 
sequence  of  covariance  estimates  having  a  straightforward 
interpretation.  If  we  form  the  matrix  Kj,**0  according  to 

K'pk  * l)  =  w;  A/*r2<*  *  X)MckWp.  (49) 

we  then  find  that 
0  •  •  0 

0  •  •  0 

0  •  •  ^**u(/*-i) 


in  a  sequence  of  estimates  of  KN.  These  have  increasing 
log  likelihoods  Z.id[K<,°>]  LJK™}  <.  •  •  •  s  Lid(K<*>] 
^  .  We  show  that  stable  points  of  this  sequence  satisfy 

the  necessary  trace  condition  for  the  maximum  likelihood 
estimate  of  Ks,  where  Kjj]  is  a  stable  point  if  Kjj+J)  -  K]P 
for  y-1,2,  ••  •. 

Each  iteration  stage  of  the  expectation-maximization 
algorithm  has  an  expectation  E  step  and  a  maximization  M 
step  that  must  be  performed  to  get  to  the  next  step.  The 
E-step  requires  evaluation  of  the  conditional  expectation 
of  the  complete  data  log  likelihood  (43)  given  the  incom¬ 
plete  data  r  and  assuming  that  the  covariance  defining  the 
complete  data  is  "Z{k). 

g[2Pl2(;']-E[Led(2,)|r,2<*'].  (44) 


where  Kj,k  +  U(i)  is  a  P  x  P  circulant  matrix  interpreted  as 
the  estimate  at  stage  k  +  1  of  the  covariance  KP(i)  of  the 
P-periodic  reflectance  process  at  delay  m  +  i.  Miller  and 
Snyder  [21]  show  that  the  (n,  m)-element  of  this  circulant 
matrix  is  given  by 

—  £  E[b(  p.i)b*((p  +  m  -  n)PJ){r.  K(Pk)\  (51) 

r  p-o 

where  (a)P  =  a  mod  P.  Equation  (51)  has  an  intuitively 
appealing  form.  If  the  reflectivity  process  b(n,i)  could  be 
observed  for  all  instants  n  =  0. 1.  •  •  • ,  P  - 1  ina  period  and 
for  each  i  independently,  then  the  maximum  likelihood 
estimate  of  the  covariance  Kp(t)  would  be  the  arithmetic 

19  - 
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average  of  the  lagged  products 
1 

-  £  b(p,i)b*((p  +  m  —  n)p,i).  (52) 

r  p- o 

Equations  (50)  and  (51)  indicate  that  one  should  simply 
substitute  the  conditional  mean  estimate  of  an  unknown 
lagged  product  into  this  expression  to  form  the  maximum 
likelihood  estimate  of  the  covariance  when  only  the  incom¬ 
plete  data  are  known. 

Estimating  2P  and  KP 

The  maximum  likelihood  estimate  of  2P  is  a  stable 
point  of  the  sequence  defined  in  (48).  The  terms  on  the 
right  side  of  this  equation  can  be  evaluated  as  follows.  Let 
the  conditional-mean  estimate  of  cP  in  terms  of  the  incom¬ 
plete  data  r  be  defined  at  stage  k  by 

#_>-E[c,|r.Zj«].  (53) 

Then,  (48)  can  be  rewritten  in  the  form 

S'**11  =  E  [(c,  -  c{k))(cP  -  c(Pk))* 1 r,  2<*>]  +  c<*>c<*)+. 

(54) 

Now  examination  of  (42)  shows  that  forming  the  condi¬ 
tional-mean  estimate  (53)  of  cP  from  r  is  a  standard 
problem  in  linear  estimation  theory.  From  Tretter  [25,  ch. 
14],  for  example,  we  find  that 

c%*>  -  2{Pk)r[r*2'Pk)r+N0i]~lr.  (55) 

Furthermore,  the  first  term  on  the  right  side  in  (54)  is  the 
covariance  of  the  estimation  error  when  cP  is  estimated 
from  r.  Also  from  Tretter  [25,  ch.  14],  we  have 

=  - 2<*)r[r+2'*)r  +  iv07] (56) 

In  summary,  the  following  steps  are  performed  to  pro¬ 
duce  a  sequence  2)?',  2)>V  •  •,  2J,*',  •  •  •  of  estimates  of  2P 
fo  which  the  corresponding  sequence  of  likelihoods  is 
nondecreasing: 

1)  set  k  =  0,  selecting  a  starting  estimate  2(/?); 

2)  calculate  the  estimate  of  cP  according  to  (55); 

3)  calculate  the  error  covariance  according  to  (56); 

4)  update  the  estimate  of  2P  according  to  (54); 

5)  if  “last  iteration”  then  stop,  else  replace  k  by  k  + 1 
and  go  to  2. 

The  starting  value  in  step  1  can  be  any  positive-definite 
diagonal  covariance  matrix  of  dimension  /*/CR  x  7*7^.  A 
sequence  of  estimates  of  Kp  having  increasing  likelihood 
is  obtained  from  the  sequence  of  estimates  of  2P  accord¬ 
ing  to  (49). 

Forming  the  Scattering-Function  Image 

From  (41),  the  diagonal  elements  of  the  7P7CR  x  /p7cr- 
dimensional  matrix  2P  are  sample  values  of  the  scattering 
function,  with  the  scattering  function  samples  at  delay 


m  +  i  given  by  the  /CR  diagonal  elements  of  the  i  th  /CR  X 
/CR-dimensional  diagonal  block  of  2P.  We  may,  therefore, 
simply  regard  2'/*  as  the  stage  k  estimate  of  the  scattering 
function.  The  stag e-k  scattering-function  image  of  the 
target  in  range  ( i  coordinate)  and  cross  range  ( p  coordi¬ 
nate)  can  be  displayed  as  follows.  Let  2 Pk)(i)  denote  the 
7 th  7cr  x  /CR-dimensional  diagonal  block  of  2(/\  and 
denote  the  /uh  diagonal  element  of  2</>(0  by  s(p,i)  = 
[o}(i)Yk)  for  p  =  0, 1  •  •  • ,  7CR-1.  Then,  s(0,  i)  is  dis¬ 
played  at  range  m  +  i  and  cross  range  corresponding  to  a 
Doppler  shift  of  zero;  j(l,/)  and  j(/CR-l,f)  are  dis¬ 
played  at  range  m  +  i  and  cross  range  corresponding  to  a 
Doppler  shift  of  o  =  l/P  and  v=—l/P,  respectively; 
s(2,i)  and  t( 7cr  - 2,  /  )  are  displayed  at  range  m  +  i  and 
cross  range  corresponding  to  a  Doppler  shift  of  v  =  2/P 
and  o  —  -2/P,  respectively;  and  so  forth,  with  s(p,  i)  and 
s( 7cr  —  p,i)  displayed  at  range  m  +  i  and  cross  range 
corresponding  to  a  Doppler  shift  of  ±  p/P  for  p  = 

•  1,2,-  •  -,(ICk  - 1)/2  when  7CR  is  odd. 

Forming  the  Reflectance-Process  Image 

It  is  interesting  to  note  that  the  k  th  stage  conditional- 
mean  estimate  of  cP,  given  the  measurements  r  and  assum¬ 
ing  that  the  second-order  statistics  of  reflectance  are  given 
by  the  k  th  stage  estimate  of  the  scattering  function,  is  used 
to  form  the  estimate  of  2P  at  stage  k  + 1  when  the 
expectation-maximization  algorithm  is  used.  This  estimate 
is  of  very  much  interest  in  its  own  right  because,  from  (36) 
and  its  definition,  the  ^CR  elements  of  cp(i)  are  the  poten¬ 
tially  nonzero  Founer-transform  coefficients  of  the  re¬ 
flectance  process  bP(i).  The  target’s  reflectance  image  at 
stage  k  is  formed  by  placing  these  elements  at  range  m  +  i 
and  cross  range  in  the  same  manner  as  described  above  for 
the  scattering- function  image. 

Convergence  Issues 

There  are  some  important  properties  of  the  iteration 
sequence  (48)  which  are  worth  mentioning.  First,  each  step 
is  in  an  improving  direction  in  the  sense  that  the  log 
likelihood  increases  at  every  step  and  continues  to  do  so 
until  a  stable  point  is  reached.  This  is  shown  by  writing 
(54)  out  as 

2</+l)  =  2y°  +  2</'©<*'2P  *  (57) 

where 

e<*>  „  r  K'k)-l{rr+  -  r+  (58) 

and  where 

/Cr('t,  =  r*2</’r  +  ;V07  (59) 

is  the  kth  estimate  of  the  covariance  K.  of  r.  Next,  the 
trace  condition  (27)  which  the  maximum  likelihood  esti¬ 
mate  must  satisfy  is  reexamined.  From  the  assumption  of 
the  P-periodicity  of  the  reflectance  process  and  the  matrix 
definitions  given,  the  admissible  variations  8K  must  be  of 
the  form 

20 


sk  -  m;  mckwpmr. 


(60) 
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Here  52  is  a  diagonal  matrix  of  the  same  dimensions  as  2. 
The  trace  condition  (27)  then  becomes 

(2Er ) 1  tr  ( K; >( rr+  -  T+ 2PT  -  N0I ) K; lT+  52  T)  -  0. 

(61) 

Using  the  fact  that  tr(AB)  =  ti(BA)  and  evaluating  this 
trace  at  the  k  th  iterate,  we  see  that  the  trace  on  the  left 
side  of  (61)  is  equal  to 

(2£r)-1tr(0(*)52).  (62) 

According  to  (57),  25>*)  is  changed  at  each  stage  by  adding 
the  diagonal  elements  of 

2<*)e(*)2<*)  (63) 

to  2</).  Define 

52<A)  -  2</)e<*)2<*)  (64) 

as  these  diagonal  elements.  Then,  evaluating  the  trace  at 
this  variation  gives 

tr(e<*>52(*>)  S:0.  (65) 

This  shows  that  the  variation  52(*’  is  in  an  improving 
direction.  Furthermore,  we  are  guaranteed  that  the  incom¬ 
plete  data  log  likelihood  is  nondecreasing  as  a  result  of  the 
M-step  of  the  expectation-maximization  algorithm  be¬ 
cause,  at  this  step,  the  conditional  expectation  of  the 
complete  data  log  likelihood,  given  the  incomplete  data 
and  the  last  iterate  for  2P,  is  maximized  over  Sp.  As 
shown  in  [20]  and  [21],  this  implies  that  the  incomplete 
data  log  likelihood  is  nondecreasing. 

Lemma  3:  Assume  that  S0  >  0  and  2(°*  is  positive  defi¬ 
nite.  Then  1)  2</)  is  positive  definite  for  all  k;  and  2)  all 
stable  points  satisfy  the  trace  condition  (27)  for  all  admis¬ 
sible  variations  (60). 

The  proof  of  the  first  part  of  Lemma  3  is  in  the 
Appendix.  For  the  second  part,  since  the  diagonal  ele¬ 
ments  of  2</)  are  positive,  (65)  holds  with  equality  if  and 
only  if  the  diagonal  elements  of  0<A)  are  zero.  Notice  that 
if  2(/  +  l)  =>  2J,*\  then  the  diagonal  elements  (64)  are  zero. 
This  implies  that  the  diagonal  elements  of  0(A)  are  zero 
and  hence  that 

tr(0<*>52)=O  (66) 

for  all  diagonal  52.  Thus  all  stable  points  satisfy  the  trace 
condition  (27)  for  admissible  variations.  From  (55),  stable 
points  of  the  sequence  2),*)  yield  stable  points  of  the 
sequence  of  conditional- mean  estimates  of  the  reflectance 
process. 

Computational  Considerations 

The  computations  required  to  produce  radar  images 
with  our  method  are  specified  by  (54)— (56).  The  number  of 
iterations  of  these  equations  required  to  produce  an  image 
near  the  convergence  point  is  presently  unknown.  Our 
experience  in  using  an  iterative  algorithm  to  produce  maxi¬ 
mum  likelihood  images  for  emission  tomography  suggests 
that  50-100  iterations  may  be  necessary,  but  this  is  only  a 


guess  that  will  not  be  verified  until  some  experiments  are 
completed.  Some  form  of  specialized  processor  to  accom¬ 
plish  each  iteration  stage  efficiently  will  probably  be  needed 
to  produce  images  in  practically  useful  times.  One  possible 
approach  is  the  following.  The  matrix  product 

T  =  {TTtMckWpMrS 

is  required  at  each  iteration  stage  and  does  not  change. 
This  /R/CR  x  N-dimensional  matrix  can  therefore  be  com¬ 
puted  once  off-line,  stored,  and  then  used  as  needed  dur¬ 
ing  on-line  computations.  Then,  at  iteration  stage  k ,  the 
following  on-line  computations  can  be  performed: 

1)  compute  the  N  X  V-dimensional  matrix  A  defined 
by/4  =  r+2;*>r  +  N0/; 

2)  compute  the  VcR  X  -dimensional  matrix  B  de¬ 
fined  by  B  =  2</*>r; 

3)  compute  BA~lr  and  the  diagonal  elements  of  2(/'  - 
BA  ~  lB  *. 

The  computations  in  3)  can  be  accomplished  in  about 
4JV  +  IRICK  -  2  time  steps  using  the  systolic  array  de¬ 
scribed  by  Comon  and  Robert  [26]  augmented,  as  they 
suggest,  by  one  row  to  accomplish  the  postmultiplication 
of  BA -1  by  r  and  by  IRICK  rows  to  accomplish  the 
postmultiplication  by  B *.  The  matrix  multiplications  in  1 
and  2  for  determining  A  and  B  can  also  be  performed 
rapidly  on  a  systolic  array.  More  study  of  implementation 
approaches  is  needed,  but  it  does  not  appear  that  the 
computational  complexity  of  our  new  imaging  algorithm 
needs  to  be  a  limitation  to  its  practical  use. 

The  choice  of  N,  IR  and  /CR  is  important  for  the 
computations.  These  parameters  are  selected  to  achieve  a 
desired  range  and  cross-range  resolution  and  are.  there¬ 
fore,  problem  dependent,  but  the  same  considerations  used 
with  other  approaches  to  radar  imaging  can  be  used  in 
selecting  them.  Choosing  the  product  lRIc R  to  be  of  the 
order  of  N  will,  in  some  sense,  make  the  imaging  problem 
well-defined  because  the  number  of  unknown  parameters 
I rIq  r  that  need  to  be  estimated  to  form  the  image  is  then 
comparable  to  the  number  of  measurements  .V.  On  the 
other  hand,  the  choice  of  P  is  unique  to  our  approach.  As 
stated,  we  need  P  ^  V,  but  no  upper  limit  is  given.  In  [24], 
it  is  shown  that  as  P  increases  toward  infinity,  so  does  the 
maximum  value  of  the  incomplete  data  log  likelihood 
function,  with  probability  one.  Thus  P  cannot  be  made 
arbitrarily  large  from  a  theoretical  standpoint.  Any  com¬ 
putation  involving  a  matrix  with  one  dimension  equal  to  P 
can  be  performed  off-line. 

IV.  Conclusion 

The  expressions  we  have  obtained  for  forming  images  of 
diffuse,  fluctuating  radar  targets  are  based  on  the  model 
stated  in  Section  II.  The  target  reflectance  is  assumed  to 
introduce  wide-sense-stationary  uncorrelated  scattering  of 
the  transmitted  signal  with  no  glint  or  specular  compo¬ 
nents  present.  The  reflectance  process  is  assumed  to  be  a 
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WSSUS  Gaussian  process  with  unknown  second-order 
statistics  given  by  a  delay-dependent  covariance  or  scatter¬ 
ing  function.  Echoes  of  the  transmitted  signal  are  received 
from  all  the  reflecting  patches  that  make  up  the  target, 
with  each  paten  introducing  some  propagation  delay, 
Doppler  shift,  and  random  amplitude-scaling  into  the  sig¬ 
nal  it  reflects.  The  superposition  of  the  echoes  from  all  the 
patches  is  received  in  additive  noise.  Thus  the  reflectance 
process  is  only  observed  indirectly,  following  a  linear  su¬ 
perposition  and  in  additive  noise;  thus  neither  the  re¬ 
flectance  process  nor  its  second-order  statistics  are  known. 
Target  images  are  made  by  displaying  estimates  of  either 
the  reflectance  process  or  its  second-order  statistics 
(scattering  function)  based  on  processing  the  received  sig¬ 
nal.  In  Section  II,  we  derived  the  trace  condition  which  the 
maximum  likelihood  estimate  of  the  covariance  of  the 
reflectance  must  satisfy,  and  we  concluded  that  this  condi¬ 
tion  is  too  complicated  to  solve  explicitly  for  the  estimate. 
This  motivated  the  introduction  in  Section  III  of  the 
incomplete-complete  data  model  and  the  use  of  the  expec¬ 
tation-maximization  algorithm,  which  results  in  a  sequence 
of  estimates  of  the  scattering  function  having  increasing 
likelihood.  A  corresponding  sequence  of  estimates  of  the 
reflectance  process  is  also  obtained. 

A  number  of  issues  have  yet  to  be  resolved  for  the 
approach  to  radar  imaging  we  have  presented.  One  of  the 
most  important  is  resolving  how  glint  and  specular  compo¬ 
nents  in  the  return  echoes  should  be  modeled  and  accom¬ 
modated  in  the  formation  of  the  images.  The  selection  of 
transmitted  signals  to  produce  good  images  is  an  impor¬ 
tant  subject  about  which  little  study  has  been  made.  The 
quality  of  target  images  obtained  with  our  new  approach  is 
not  known  at  present;  to  study  this  issue,  we  are  presently 
implementing  a  computer  simulation  so  that  comparisons 
to  alternative  processing  strategies  cam  be  made.  The  equa¬ 
tions  we  have  developed  are  computationally  demanding 
because  the  matrices  involved  can  be  of  large  dimension 
and  the  iteration  must  be  performed  repeatedly  until  a 
stable  point  is  approached.  It  is  therefore  important  to 
determine  the  conditions  under  which  our  approach  yields 
radar  images  of  sufficiently  improved  quality  compared  to 
existing  approaches  to  warrant  the  development  of  special 
processing  architectures  that  will  make  it  practical.  The 
computer  simulations  should  be  of  some  help  in  this.  At 
this  time,  however,  the  only  experimental  results  suggest¬ 
ing  the  efficacy  of  our  method  are  those  reported  in  [21  j 
and  [22]  for  estimating  power-density  spectra  in  one  di¬ 
mension. 
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Appendix 

Proof  of  Lemma  1 

From  the  definition  of  the  log  likelihood  function  in  (26),  we 
have 

+  a&K)~  L{i(K)) 

(X 

-  -  i  r*  ((  Kr  +  a2ETS*  SfCS )  '  1  -  Kf 1 )  r 

-  ~  ln(det(  Kr  +  a2£>S*  SfCS)  det(  fCf1))  (Al) 

where  Kr  is  the  covariance  of  the  incomplete  data  r  as  given  in 
(26).  Examining  the  first  term  on  the  right,  we  have 

-ir+ATr-1((/  +  a2Fr5+  SfCS  Kfl)~l  -  i)  r 

-  ^  r*  K; 1  (  al  ErS*  SKSfCf  l)(/  +  a2  ErS~  SfCSfC;l)~'r 

-r*K;l2ETS*  SfCS K;'r  + 0(a).  (A 2) 

Examining  the  second  term  on  the  right  in  (Al),  we  have 

—  ^  ln(det(  /  +  a2ErS~  SKSfCf1)) 

1 

- - ln(det(  /  +  aB)) 

a 

i 

” - ln(l  +  atr(fl)  +  +o"det(S)) 

a 

-~t r(B)  +  0(a)  (A3) 

where  B  is  defined  in  the  first  equality.  For  any  fC  6  Q  to  be  a 
local  maximum,  a  small  variation  in  AT  in  an  admissible  direction 
cannot  increase  Li(i(K),  or 

1 

lim  -(L^K  +  aSK)- Lid(K))  <0  (A4) 

aiO*  a 

for  all  admissible  variations  SfC.  If  AT  is  an  interior  point,  then 
-  SfC  is  also  an  admissible  variation  and  (A4)  becomes  an 
equality.  Substituting  (A2)  and  (A3),  we  get 

rK;'2ETS*  SfCS K;'r-u(2ETS~  SfCS K;1)  -0.  (A5) 

which  is  the  trace  condition  (27). 

Proof  of  Lemma  2 

Suppose  that  K  satisfies  (27)  and  (29)  for  all  admissible 
variations  SfC.  We  now  show  that  (29)  is  simply  the  second 
derivative  of  Z.id(Af)  in  the  direction  SfC  by  taking  the  limit 

lim  —  u((2ErS*(  K  +  a8fC)S+  N0I)~' 

<*  —  0  [ot( 

■(rr*  -  N0I  -2ETS*  (  K  +  aSK) S) 

(2 EtS~KS  +  a2ErS*  SfCS  +  ,V0/)  ~  '  S*  SfCS 

-  fC;l(  rr *  -  ff0I  -  2£rS*  KS)  Kf  lS~  SKS) 

-  ti(K;'2ErS~  8fCSK;'(2ErS*  KS 

+  /V0/-2rr”  )  Kfl2ErS*  SKS).  (A6) 
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Thus  the  conditions  of  Lemma  2  are  the  standard  sufficient 
conditions  for  a  point  K  to  be  the  local  maximum  of  Lld(K). 
Equation  (27)  says  that  the  first  derivative  of  Lid(K)  is  zero  at 
K.  Equation  (29)  says  that  the  second  derivative  is  negative 
definite  along  admissible  variations  from  K.  A  necessary  condi¬ 
tion  for  K  to  be  a  relative  maximum  is  that  this  last  expression 
evaluated  at  K  to  be  equal  to  or  less  than  zero  for  all  admissible 
variations  SK.  Under  the  assumptions  in  Section  IV,  admissible 
variations  are  given  by  (60).  Substituting  (60)  into  (A6)  and 
evaluating  for  all  diagonal  matrices  <52  gives  the  second-order 
necessary  condition. 

Estimating  a  General  Toeplitz  Matrix 

In  Section  IV,  we  derived  a  sequence  of  estimates  for  a 
covariance  matrix  subject  to  the  constraint  that  the  estimates  be 
circulant  Toeplitz  matrices.  For  completeness,  we  develop  and 
discuss  the  equations  for  estimating  a  covariance  matrix  subject 
to  the  weaker  constraint  that  the  estimates  be  general  Toeplitz 
matrices.  Similar  equations  for  other  constraints  on  the  Toeplitz 
matrices  are  easily  obtained  by  mimicking  the  steps  in  the  main 
body  of  this  paper. 

Let  the  complete  data  be  {b,w},  and  let  b  be  normally 
distributed  with  zero  mean  and  covariance  K,  as  given  in  (27). 
The  complete  data  log  likelihood  is 

£*,(*»  --la(det(Ar,))-6;^'^  (A7) 

where  all  terms  that  are  not  a  function  of  KP  have  been  sup¬ 
pressed.  Maximizing  this  function  gives  the  trace  condition 

u{K-'{bb*-K)K-l8K)~Q,  (A8) 

which  the  maximum  likelihood  estimate  K  must  satisfy.  Perform¬ 
ing  the  E-  and  M-steps  of  the  expectation-maximization  algo¬ 
rithm  yields  the  following  iteration  sequence  for  the  elements 
K(n,i),  n  -  0,1,- •  •,  Af-1,  of  the  covariance  matrix  K(i)  de¬ 
fined  in  (23): 

i  r  v  -  n  - 1 

- E  21  b(j,m  +  i)b*(j  +  n,m  +  i)\r,K{k)  . 

N~n  i- o 

(A9) 

In  matrix  form, 

A-'**1'  -  K{k)  +2ETK(k)SK<,k)~l 

(rr* -lErS*K{k)S- N0l)  K'k'-'S~  K(k)  (A10) 

where 

fC*k)  -  2ETS* Kik)S  +  N0I.  (All) 

If  this  iteration  converges  to  a  stable  point,  then  the  trace 
condition  is  satisfied  at  this  point,  as  may  be  shown  by  using  the 
same  arguments  as  in  Section  IV.  It  is  worth  restating  that  the 
reason  this  iteration  is  not  recommended  here  is  that  the  proba¬ 
bility  that  the  iteration  sequence  generates  a  singular  estimate  for 
K  approaches  one  as  <V  gets  large.  By  restricting  consideration  to 
Toeplitz  matrices  with  circulant  extensions,  the  log  likelihood 
function  is  bounded  with  probability  one  for  finite  extensions 
and  a  positive  definite  K  is  generated  with  probability  one,  as 
proven  by  Fuhrmann  and  Miller  (24). 


Proof  of  Lemma  3,  Par t  1 

Assume  that  the  initial  guess  2)?’  for  is  positive  definite 
and  that  V0  >  0.  We  will  now  show  that,  if  2),*1  is  positive 
definite,  then  so  is  2p**l\  and  thus,  by  induction,  2)>* 1  is 
positive  definite  for  all  k.  One  key  to  following  this  derivation  is 
the  matrix  identity 

B(I  +  AB)~l -(/+  BAy'B.  ( A12) 

This  identity  is  used  to  rewrite  (57)  as 

2<*  +  ‘>  -  //(  V02(/)rr*  2J,*>  +  V022J,*> 

+  2{Pk)Trr*r*VPk))H* 

-n0(h  2,/»r)(r*2<,*,/ir*) 

+(ff2V*,i7)(r*r*2,/,/r ) 

+  V02//2',‘)7/*  (A13) 

where  we  have  defined  H  according  to 

ff-(2'/.*Tr*V0/)'1.  (A14) 

Dearly,  all  the  diagonal  elements  of  (A13)  are  greater  than  or 

equal  to  zero.  To  show  that  they  are  strictly  positive,  we  look  at 
the  last  term  and  get  that  the  i  th  diagonal  element  is 

y-o 

“  ‘Vo  *  £  (A15) 

y-o 

which  is  clearly  positive  when  ,V0  >  0  since  H  is  invertible  and  all 
diagonal  elements  of  2J,* ’  are  positive. 
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ABSTRACT 

We  have  developed  a  new  approach  to  maximum-likelihood  spectrum  estimation 
of  wide-sense  stationary  processes  from  noisy  data.  A  statistical  model  for  the  data  is 
defined.  The  process  whose  spectrum  is  sought  is  wide-sense  stationary,  periodic  and 
Gaussian,  and  its  observations  are  corrupted  by  an  additive  white  noise.  A  maximum- 
likelihood  formulation  of  this  problem  has  been  derived,  and  the  equations  are  solved 
numerically  via  the  expectation-maximization  algorithm.  This  approach  presents  several 
attractive  features,  an  important  one  being  that  the  noise  corrupting  the  observations  is 
now  taken  into  account. 

We  present  some  recent  developments  for  this  problem.  The  statistical  perfor¬ 
mance  of  the  new  maximum-likelihood  spectrum  estimator  is  studied  both  theoretically 
and  numerically.  Comparison  with  traditional  estimators  such  as  the  periodogram 
highlight  several  strong  points  of  the  method.  We  also  identify  certain  limitations, 
namely  the  instability  of  estimates  for  high  noise  levels.  These  limitations  can  be  allevi¬ 
ated  if  a  priori  information  about  the  signal  is  available.  Two  such  problems  are  dis¬ 
cussed  in  which  the  information  at  hand  has  the  form  of  a  constraint  on  the  input  signal- 
to-noise  ratio. 

We  show  how  such  information  can  be  incorporated  in  the  maximum-likelihood 
estimation  procedure.  First  we  assume  the  signal  power  to  be  known.  Theoretical  issues 
of  existence  and  uniqueness  of  the  solution  are  discussed.  We  proceed  with  a  problem  in 
which  the  information  is  less  complete,  when  only  an  upper  bound  and/or  a  lower  bound 
on  the  signal  power  are  available.  The  statistical  performance  of  both  constrained  esti¬ 
mators  is  quantitatively  studied. 


*  This  work  was  supported  by  contract  number  N00014-86-K-0370  from  the  Office  of  Nival  Research. 
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1.  Introduction 

A  promising  approach  to  maximum-likelihood  estimation  of  Toepli'7.  constrained  covariance 
matrices  has  been  proposed  recently  [1].  Several  further  developments  cai  je  considered.  First,  this 
method  also  applies  to  the  dual  problem  of  spectrum  estimation.  Another  issue  of  interest  is  that  the  statist¬ 
ical  model  can  account  for  the  presence  of  additive  noise  corrupting  the  observations  and  for  linear 
transformations  of  the  process  whose  covariance  or  spectrum  is  sought.  These  considerations  have 
motivated  a  new  approach  to  high-resolution  delay-doppler  radar  imaging,  where  a  major  goal  is  to  pro¬ 
duce  estimates  of  the  target’s  scattering  function  [2],  In  the  special  case  of  a  point  target  and  a  constant 
envelope  transmitted  signal,  this  reduces  to  a  spectrum  estimation  problem. 

This  paper  describes  some  recent  developments  for  this  problem.  We  study  the  statistical  perfor¬ 
mance  of  the  new  maximum-likelihood  spectrum  estimator  both  theoretically  and  numerically.  Com¬ 
parison  with  traditional  estimators  such  as  the  periodogram  highlight  several  strong  points  of  the  method. 
We  also  identify  certain  limitations,  namely  the  instability  of  estimates  for  high  noise  levels.  These  limita¬ 
tions  can  be  alleviated  if  a  priori  information  about  the  signal  is  available.  Two  such  problems  are  dis¬ 
cussed  here  in  which  the  information  at  hand  has  the  form  of  a  constraint  on  the  input  signal-to-noise  ratio. 

This  paper  is  organized  as  follows.  Our  model  is  presented  in  Section  2.  A  maximum-likelihood 
formulation  of  the  problem  is  given  in  Section  3,  and  the  equations  are  solved  via  the  expectation- 
maximization  algorithm.  Section  4  is  devoted  to  a  statistical  performance  analysis  of  this  estimator  and  a 
comparison  with  two  other  methods.  In  Section  5  we  show  how  a  priori  information  on  the  signal  can  be 
incorporated  in  the  maximum-likelihood  estimation  procedure.  First  we  assume  the  signal  power  to  be 
known.  Theoretical  issues  of  existence  and  uniqueness  of  the  solution  are  discussed.  Section  5  deals  with 
a  less  complete  knowledge,  where  only  an  upper  bound  and/or  an  lower  bound  on  the  signal  power  are 
available.  The  last  section  is  devoted  to  a  quantitative  study  of  the  statistical  performance  of  both  con¬ 
strained  estimators. 

2.  Mode! 

The  following  is  derived  from  the  model  presented  in  [1]  for  a  point  target  and  a  constant  envelope 
transmitted  signal.  The  observation  is  an  N-vector  sample  of  a  wide-sense  stationary,  periodic,  Gaussian 
process  corrupted  by  an  additive  noise : 

r  =  b  +  w  ,  (2.1) 

where  b  contains  N  consecutive  samples  of  a  zero-mean  periodic  process  bp  with  length  P  >  N,  and  w  is 
an  zero-mean  white  Gaussian  noise  with  variance  N0,  uncorrelated  with  b .  The  periodicity  assumption  is 
required  to  guarantee  that  the  likelihood  function  is  bounded  above;  therefore,  there  exists  the  maximum- 
likelihood  estimator  [1], 

•  This  work  was  supported  by  contract  number  N00014-86-K-0370  from  the  Office  of  Naval  Research. 
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Now  we  define  the  spectral  process  associated  with  bp  to  be  the  DFT  of  one  period  of  bp.  Assume 
that  we  are  interested  in  estimating  only  M  of  the  components  of  this  spectral  P-vector  (1  <  M  <  P),  the 
other  components  being  zero  with  probability  1;  let  c  be  this  M-vector.  This  assumption  is  introduced  to 
deal  with  the  bandlimited  spectra  encountered  in  radar  applications,  which  arise  because  radar  targets  have 
finite  extent  [2].  c  is  a  Gaussian  random  M-vector  with  diagonal  covariance  Z,  whose  entries  c r(i),  i  = 
are  real  and  positive,  c  and  b  are  related  by  a  linear  transformation  : 

b  =  T*c  ,  (2.2) 

where  we  have  defined  the  M  x/V  matrix  T,  consisting  of  the  first  N  rows  and  the  outer  M  columns  of  the 
PxP  DFT  matrix.  The  superscript  t  denotes  the  Hermitian-transpose  operator  on  matrices.  Our  model  for 


the  observations  can  now  be  written  as 

r  -  r'c  +  w  . 

(2.3) 

The  covariance  matrix  for  r  is  given  by 

Kr  =E [rrf]  =  TfZT  +  jVq/y  , 

(2.4) 

where  /iV  is  the  jVxjV  identity  matrix. 

3.  Spectrum  Estimators 

In  this  section  we  introduce  a  maximum-likelihood  spectrum  estimator  for  the  model  (2.3),  denoted 
by  ML1.  We  also  define  two  estimators  which  will  be  analyzed  and  compared  to  ours  in  the  next  section. 
The  first  one  is  the  maximum-likelihood  estimator  derived  assuming  noise-free  data,  denoted  by  MLO  ;  the 
second  one  is  the  periodogram. 


3.1.  ML1  Estimator 

From  (2  l),  the  likelihood  function  for  Z  is 

L  (r  ,Z)  =  -‘/z  In  det  (T^ZT  +  NqIn)  -  Vz  rf(T^Zr  +  NcJsl)~lr  .  (3.1) 

Maximizing  the  likelihood  with  respect  to  Z  yields  the  necessary  trace  condition  which  the  estimate  Z  must 
satisfy  [1(2]: 

Tr  [r(Tftr  +  NoTvr'irr'-r'ir  -  Ay,v)(rfzr  +  AV.vr’r'SZ]  =  0  ,  (3.2) 

for  all  M xAf  diagonal  matrices  8Z.  This  trace  condition  is  a  nonlinear  equation  in  Z.  Generally  it  cannot  be 
solved  directly  in  closed-form,  so  some  numerical  search  procedure  must  be  implemented.  An  elegant  solu¬ 
tion  is  the  expectation-maximization  (EM)  algorithm  used  in  [1,2].  An  initial  estimate  H0)  is  selected.  At 
step  k+1  (k  =  0,1,..)  the  estimate  is  updated  according  to 

Z(*+1)  =  argmax  Q  (Z  j  t'*5)  (3.3) 

where 


Q  (Z|Z(i>)  =  -V:  *£  h  <**(0  ~  S' 


i=0 


i= 0 


(3.4) 


and 


e  [  )c  (i)j2|  r  j5l)]  =  [f*’  -  z(*)nrrza,r+Af(/w)-1r^z(*>  +  z(*)r(T'Z(t)r+/V(/iV)-1 

xrrf(r^Z(i)r+/V0//vr1rff*)]i,-  .  (3.5) 


This  algorithm  produces  a  sequence  of  estimates 


<r(/)<t+1)  =  £[|c(i)|2|r,Z(*)] 


(3.6) 


having  increasing  likelihood.  It  can  be  shown  that  the  stable  points  of  this  algorithm  satisfy  the  necessary 
trace  condition  for  a  maximizer  [2].  The  issue  of  uniqueness  is  addressed  in  [3], 
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Special  case  :  N  =  M  =  P 

In  this  special  case,  a  full  period  of  the  process  is  observed,  r  is  equal  to  the  P  xP  DFT  matrix  \Vp . 
A  closed-form  expression  for  t  can  be  derived: 

<?(D  =  max(0,|(Wpr)(O|2-  N0)  .  (3.7) 


3.2.  MLO  Estimator 

Additive  noise  corrupting  observations  is  usually  not  included  separately  in  approaches  to  spectrum 
estimation.  This  model  was  assumed  in  [1].  The  sequence  of  estimates  of  I  is  still  given  by  (3.6)  and  (3.5), 
in  which  we  now  let  N0  =  0.  We  call  this  the  MLO  estimation.  Clearly  MLO  and  ML1  are  equivalent  for 
noise-free  problems. 

Special  case  :  N  =  M 

The  problem  for  which  the  number  of  observations  (N)  is  equal  to  the  number  of  parameters  to  be 
estimated  (M)  is  of  some  practical  interest.  It  also  turns  out  that  the  trace  condition  can  be  solved  in 
closed-form  in  this  instance.  The  matrix  f  is  then  invertible,  indicating  the  existence  of  a  one-to-one  map¬ 
ping  between  r  and  c.  The  MLO  estimator  is  simply 

a2(/)=|(rV)(i)|2,  (3.8) 

where  I~f  denotes  (T-1/. 

3.3.  Periodogram 

The  periodogram  estimate  of  the  spectrum  is  defined  as  the  (scaled)  magnitude-squared  Fourier 
transform  of  the  N  observations  padded  with  P-N  zeroes  [3].  The  first  M  spectrum  samples  are  then  given 
by 

<?(0  =  (PW)|(Tr)(0|2  (3.9) 

Special  case  :  N  =  A/  -P 

When  N  =  M  =  P,  the  periodogram  and  MLO  estimates  are  the  same. 

4.  Bias  Performance  Analysis 
4.1.  Performance  Evaluation 

In  this  section,  we  estimate  Z  for  the  model  (2.3)  and  study  the  statistical  performance  of  the  three 
estimators  above.  For  each  method  the  bias  is  evaluated,  where 

Bias  [Z]  =  E  [Z]  -  Z  (4.1) 

As  we  shall  see  in  Section  4.3,  the  performance  strongly  depends  upon  the  input  signal  to  noise  ratio 
defined  by 

SNRUl=E0/N0,  (4.2) 

where  E0  is  the  average  power  of  the  process,  defined  by 

£0  =  (1/P)7>  [Z],  (4.3) 

In  sections  4.2  and  4.3,  we  evaluate  the  bias  for  the  estimators  derived  in  Section  3.  Whenever 
closed-form  expressions  are  not  available,  computer  simulations  are  performed.  Typically  3000  realiza¬ 
tions  are  generated  for  each  process.  For  a  given  estimator,  (4.1)  is  then  estimated  from  the  3000  esti¬ 
mates.  The  analysis  is  carried  out  at  various  input  SNR  levels.  Much  effort  was  made  for  the  special  case 
M  =  N,  This  provides  insight  into  the  problem  since  the  MLO  equations  can  be  solved  in  closed  form.  The 
choice  of  P  is  free,  so  long  as  P  5  N  [2]. 
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4.2.  Closed-form  Expressions  for  Estimator  Bias  Performance 

(a)  MU 

As  indicated  in  Section  3.1,  no  closed-form  expression  for  the  estimator  is  available,  so  the  evalua¬ 
tion  of  the  bias  is  obtained  by  computer  simulation. 

(b)  MLO 

Closed-form  expressions  for  MLO  can  be  derived  when  M  =  N.  The  results  are  presented  below. 


Combining  (2.3)  and  (3.8),  we  can  write 

&(i)=\(c  +rtw)(i)\2. 

(4.4) 

Taking  the  expectation  of  (4.4),  we  get 

£[a2(i)]  =  a2(i)  +  N0(TT'r1,  , 

(4.5) 

which  implies 

Bias[cr1(i)]=NQ(TTtTXu  . 

(4.6) 

The  bias  is  due  to  the  noise  corrupting  the  observations  and  is  proportional  to  its  variance.  The  sensitivity 
of  the  bias  to  the  noise  is  determined  by  the  diagonal  entries  of  the  matrix  (TT')'1 . 


(c)  Periodogram 

Bias 

Combining  (2.3)  and  (3.9),  we  write  the  periodogram  estimates  in  the  equivalent  form 

o*(Q  =  (P  in)  icrr'c  +  rw)(o|2 .  (4.7) 

Taking  the  expectation  of  (4.7),  we  get 

£[<?(<■)]  =  (.PIN)  (Tr'XTT^  +  N0  ,  (4.8) 

and 

Bias  [<?({)]  =  ( (PIN)  (TT'ZIT'X  -  <T(Q  ]+N0.  (4.9) 

The  bias  contains  two  terms.  The  second  is  due  to  the  noise  and  is  proportional  to  N0.  The  other  term  is 
independent  of  N0.  Even  for  noise-free  observations,  the  periodogram  is  a  biased  estimator  of  E  unless  IT' 
is  the  identity  matrix.  This  would  be  the  case  only  for  N  =  M  =  P  (observation  of  a  full  period  of  the  pro¬ 
cess)  or  N/M  — ♦  oo  (infinite  data). 

4 3.  Simulation  results 
Process  1 

The  first  process  we  consider  is  real  and  has  period  P  =  10.  Its  spectrum  is  symmetric  and  lowpass  (M  =  5). 
All  nonzero  spectrum  samples  are  identical : 

o2(i)=  1,  i  =  0,..,4  . 

The  number  of  observations  is  N  =  M  =  5.  The  noise  variance  N0  ranges  from  0  to  1.  Figure  1  shows  the 
bias  for  the  estimators  of  <T(2)  as  a  function  of  SNR^,  according  to  the  definitions  (4.1)  and  (4.2).  In  the 
absence  of  additive  noise  (  SNRin  — »  «»),  ML1  and  MLO  are  the  same.  Both  are  unbiased  estimators.  The 
periodogram,  however,  is  biased.  When  N0  increases  from  0,  the  performance  of  the  estimators  is  roughly 
constant  so  long  as  SNRin  remains  above  some  threshold.  For  larger  N0,  all  three  estimators  exhibit  a  strong 
degradation  in  performance.  Comparing  the  thresholds  for  MLO  and  ML1,  we  see  the  tremendous  improve¬ 
ments  brought  by  taking  the  noise  into  account  in  the  model.  Typically,  for  a  same  bias,  ML1  will  have  the 
same  performance  as  MLO  operating  in  a  20  dB  noisier  environment. 
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Process  2 

As  shown  in  the  next  section,  the  periodogram  does  not  perform  well  for  rough  spectra.  This 
motivated  our  study  of  a  sharply  peaked  spectrum.  The  process  has  period  P  =  10,  and  a  single  nonzero 
spectrum  component 

<r(0)=I. 


There  is  justN  =  M  =  1  observation. 

The  bias  for  the  estimators  of  cr(0)  is  plotted  as  a  function  of  SNRin  in  Figure  2.  In  the  absence  of 
additive  noise,  the  bias  of  the  periodogram  is  -90%  of  o^O).  Clearly,  the  conventional  estimator  is  outper¬ 
formed  by  MLO  and  ML1.  It  should  also  be  noticed  that  for  this  process,  the  improvement  of  ML1  over 
MLO  is  quite  reduced. 

Computational  Considerations 

The  convergence  rate  of  the  EM  algorithm  depends  on  several  parameters.  The  computation  time  for 
each  iteration  is  of  order  Af  N2.  The  number  of  iterations  required  for  convergence  of  the  algorithm  grows 
as  M  and  N  increase.  For  ML1,  more  iterations  are  required  as  N0  increases,  especially  in  the  threshold 
region  and  beyond.  Typical  figures  are:  for  process  1  with  N0  =  0.1,  30  iterations  are  required  before  the 
spectrum  estimates  are  stable;  when  jV0=  1,  300  iterations  must  be  performed.  Our  algorithm  is  imple¬ 
mented  on  a  Masscomp  model  5500.  Running  the  program  on  3000  realizations  in  the  latter  case  is  typi¬ 
cally  completed  in  6  CPU  hours.  We  are  presently  implementing  these  algorithms  on  a  mesh-connected 
1024  processor  (DAP  by  Active  Memory  Technology),  and  we  expect  a  major  reduction  in  the  time 
required  to  produce  estimates. 


4.4.  Discussion 

The  results  derived  above  suggest  additional  comments  on  the  periodogram.  It  can  be  shown  that 
(4.8)  can  be  written  in  the  alternative  form  [4] 

E [<?(/)]  =  ofy)  *  aN(i)  +  N0 ,  (4.10) 

where  *  denotes  the  discrete  convolution  operation,  and  £2w(z)  is  the  DFT  of  the  window 

:|n|<A/  for?  >2N  (4.11) 


cow(«)=l-^ 


COjv(n) 


=  0 

N 

_  P_ 

~  N 


for  2N  >  P  >  N 


: P-N  <|n  [<  y 


The  main  lobe  of  £2*  has  width  — .  Consider  now  a  process  made  of  sharp  isolated  spectral  peaks,  such  as 

N 

Process  2  above.  Equation  (4.10)  shows  that  the  energy  in  these  peaks  is  spread  out  as  a  result  of  the  con- 

p 

volution  operation  ;  consequently,  these  peaks  are  grossly  estimated  when  —  is  large. 

N 


5.  Constrained  maximum-likelihood  estimation 
5.1.  Description  of  the  problem 

In  this  section  we  focus  our  attention  on  ML1.  An  examination  of  Figures  M  suggests  that  ML1 
suffers  in  certain  situations.  When  SNR ^  is  low,  the  estimates  are  biased  ;  as  we  shall  soon  sec,  their  vari¬ 
ance  is  also  large.  Although  the  maximum-likelihood  estimator  is  asymptotically  unbiased  and  efficient, 
these  properties  are  not  guaranteed  in  the  small-sample  problems  considered  in  Section  4.  This  limitation 
can  be  alleviated  if  a  priori  knowledge,  such  as  SNRin,  is  available.  Since  N0  is  known,  such  a  constraint 
on  the  signai-to-noise  ratio  can  be  translated  into  a  constraint  on  the  signal  power  that  must  be  sausficd  by 
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the  maximum-likelihood  estimates.  Now  we  show  how  this  constraint  can  be  incorporated  into  the  EM 
algorithm.  The  constrained  estimates  exist  and  are  unique. 

In  Section  5.2,  SNRin  is  known.  In  Section  5.3,  our  knowledge  is  more  incomplete,  and  only  an  upper 
bound  and/or  a  lower  bound  on  SNR u  are  available. 

5.2.  Known  SNRM 

The  equations  for  ML1  presented  in  Section  3.1  can  be  modified  as  follows  to  satisfy  the  constraint 
At  each  step  of  the  EM  algorithm,  we  maximize  Q  (2|  2(t>)  defined  in  (3.4),  subject  to  the  power  constraint 

£l(?(Q  =  P£0  =  S  .  (5.1) 

;=o 

where  E0  is  the  signal  power.  The  solution  also  maximizes 

2(Z|£VMIIo2(0-S),  (5.2) 

;=o 

where  X.  is  a  Lagrange  multiplier.  Taking  the  gradient  of  (5.2)  with  respect  to  2,  we  obtain  a  quadratic 
equation  for  each  spectral  component 

2  X  c4(i)  -  o2^)  +  C,  =  0  ,  (5.3) 

where 

C,=£[|c(i)|2|r  v(‘)] 

is  calculated  according  to  (3.5).  The  solution  to  (5.3)  is 


a2(i)  = 


1  +  /,Vl-8C,X 


:X*0 


:X  =  0, 


where  /,  is  either +1  or  -1.  The  equation  for  X  is 


m-i  - - 

4SX-M  =  2 /.V1-8C.X  .  (5.5) 

i=0 

In  general  this  nonlinear  equation  in  X  cannot  be  solved  in  closed-form.  Furthermore,  an  ambiguity  subsists 
about  the  choice  of  the  signs  /,• .  The  latter  problem  is  solved  by  application  of  the  following  theorem  : 

Theorem 

Assume  that  C0  >  C( ,  i  =  1...JM-1.  Then 
(1) 

/,  =  -1  = 

/0  =  + 1  :  S  <  2C0  [M  -  I’/.Wl^C./Co)] 

i=i 

=  -1  :  else 

(2)  X  is  the  largest  nonzero  solution  of 

(4 SX-M  +  I'/.WHC, /Cq))2=1-8Co^  for  5  ,  (5.6a) 

1*1 


X  =  0,  for  5=2  Cj . 

is! 
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X  is  upper-bounded  by  1/8C0,  and  (5.6a)  can  be  solved  numerically  for  X.  Note  that  the  particular  case 
(5.6b)  is  also  the  solution  to  the  unconstrained  maximization  problem.  Next,  a 2(i)(i+1)  is  calculated  from 
(5.4).  The  whole  procedure  is  repeated  at  each  maximization  step  of  the  EM  algorithm.  Note  that  because 
of  the  highly  nonlinear  nature  of  the  problem,  no  analytic  expression  is  available  for  the  constrained  esti¬ 
mator,  even  in  the  special  case  mentioned  in  (3.7). 

5_3.  Known  upper/lower  bound  on  SNRin 

In  this  section,  the  a  priori  knowledge  about  SNRin  has  the  form  of  an  upper  bound.  Our  approach 
parallels  that  of  the  previous  section,  with  the  upper  bound  now  expressed  as  an  inequality  constraint  on  the 
estimated  signal  power.  At  each  step  of  the  EM  algorithm,  we  maximize  Q  (Z|Z(i))  defined  in  (3.4),  sub¬ 
ject  to  the  inequality  constraint 

M£?(i)<PEmtx  =  Sm„,  (5.7) 

/=o 

where  £mix  is  the  upper  bound  on  the  signal  power.  If  the  unconstrained  solution  satisfies  the  upper  bound, 
the  constraint  is  inactive  and  the  estimate  is  given  by  (3.6).  Otherwise,  the  constraint  is  active,  and  as  in 
Section  5.2,  the  solution  is  the  maximizer  of  the  expression  (5.2). 

We  can  expect  the  performance  of  this  estimator  to  be  strongly  conditioned  by  the  choice  of  £mlv  In 
the  limiting  case  £mix  — >  <»,  the  constraint  is  always  inactive  and  the  estimator  is  equivalent  to  the  uncon¬ 
strained  estimator.  For  the  other  extreme  case  £mJI  -»  0,  the  constraint  is  always  active.  A  lower  bound  or 
simultaneous  upper  and  lower  bounds  are  treated  in  exactly  the  same  manner. 

6.  Simulation  results 

In  this  section,  we  apply  the  SNR-constrained  estimators  derived  above  to  Process  1,  and  we  evalu¬ 
ate  numerically  both  their  bias  and  mean-squared  error,  where 

Var  [Z]  =  £  [Z2]  -  (£  [Z])2  (6.1) 

USE  [Z]  =  £  [(Z-Z)2]  =  Var  [Z]  +  ( Bias  [Z])2  .  (6.2) 

The  output  signal  to  noise  ratio  matrix  is  defined  as  follows  : 

SNR0U1  [t ]  =  E  [X]  (USE  [Z])-*  .  (6.3) 

Figures  3  and  4  give  a  plot  of  the  bias  and  SNRoia  for  three  different  estimators  of  <r(2)  as  a  function 
of  SNRin,  according  to  the  definitions  (4.1),  (4.3),  and  (6.3).  The  estimators  represented  on  these  figures 
are:  the  two  constrained  estimators  of  Section  5,  respectively  denoted  by  EQ-MLE  and  INEQ-MLE,  and 
defined  for  the  power  constraint  5=5  and  5,^,  =  15,  respectively  ;  and  the  unconstrained  estimator  ML1 
of  Section  3.1. 

In  the  absence  of  additive  noise  (  SNR^  -»  »),  ML  I  and  EQ-MLE  are  unbiased.  However,  the 
periodogram  is  biased,  and  so  is  INEQ_MLE.  For  the  latter,  this  can  be  understood  as  follows.  The  distri¬ 
bution  of  the  sum  S  of  the  M  estimates  ofy)  is  truncated  (to  5™*  =  15),  and  therefore  the  sum  of  all  biases 
is  negative. 

When  N0  increases  from  0,  the  performance  of  the  estimators  is  roughly  constant  so  long  as  SNRm 
remains  above  some  threshold.  For  larger  N0,  all  estimators  exhibit  a  degradation  in  performance.  Note 
that  for  the  SNR-constrained  estimators,  each  bias  is  upper-bounded  by  Sc  -  and  lower-bounded  by 
-  ofy),  where  Se  is  the  constraint.  Comparing  the  SNRoul  performance  in  Figure  2,  we  see  other  favor¬ 
able  effects  of  incorporating  SNR  constraints  into  the  problem.  For  low  N0,  SNRout  is  improved.  This  is 
due  to  the  estimates  having  a  lower  variance,  which  is  the  dominant  term  in  SNR0UJ.  For  very  noisy  data, 
the  performance  of  the  estimators  is  clearly  improved.  We  can  easily  derive  a  lower  bound  for 
SNR^ird)] : 

- - __<i. 

max  [Sc  -<r(i)  ,  cr(i)] 
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This  bound  is  independent  of  No- 
Conclusions 

In  this  paper,  we  have  described  an  approach  to  spectrum  estimation  from  noisy  data,  based  upon  a 
statistical  model  for  the  observations.  First  we  derive  a  maximum-likelihood  estimator,  and  evaluate  its 
statistical  performance.  A  comparison  is  made  with  two  other  methods  that  do  not  take  the  additive  noise 
into  account  One  is  the  traditional  periodogram  and  the  other  is  the  maximum-likelihood  estimator  derived 
for  a  noise-free  model.  It  is  shown  that  the  new  estimator  offers  a  much  better  bias  performance.  The 
improvement  over  the  periodogram  is  particularly  noticeable  for  rough  spectra:  The  bias  of  the  periodo¬ 
gram  was  as  high  as  90%  for  the  process  #2  we  considered. 

In  general  however,  the  maximum-likelihood  estimates  are  still  unstable  at  high  noise  levels.  In  the 
second  step  of  our  study,  we  refine  our  technique  to  improve  the  performance  when  some  side  information 
exists.  We  have  studied  one  such  problem  in  which  some  information  about  the  signal-to-noise  ratio  is 
available.  The  performance  for  the  SNR-constrained  estimators  has  been  numerically  evaluated,  and  com¬ 
pared  with  that  of  the  unconstrained  estimator  and  of  the  periodogram.  The  new  estimators  perform 
significantly  better  than  their  competitors  for  low  SNRin.  Because  of  the  SNR  constraint,  the  estimates  are 
not  allowed  to  take  on  the  large  values  that  were  produced  in  the  unconstrained  estimation  problem.  This 
results  in  the  estimates  having  a  lower  variance.  One  additional  feature  of  our  approach,  and  an  attractive 
one,  is  its  versatility.  Only  a  slight  modification  of  the  (unconstrained)  algorithm  is  required. 
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Figure  3.  Bias  (d2(2)) 

ML1  (solid  line),  EQ-MLE  (5=5,  dotted  line),  INEQ-MLE  (5mlx=l5,  dashed  line) 
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Figure  4.  SNR0Ut  (cT(2)) 

ML1  (solid  line),  EQ-MLE  (5=5.  dotted  line),  INEQ-MLE  (5m,x=l5,  dashed  line) 
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6.3  Reprint  of:  J.  A.  O’Sullivan,  D.  L.  Snyder,  and  P.  Moulin:  "The  Role  of  Spectrum  Esti¬ 
mation  in  Forming  High-Resolution  Radar  Images",  Proc.  ICASSP  1989,  Glasgow,  U.K., 
May  1989. 
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ABSTRACT 

We  have  developed  a  new  approach  to  forming  high-resolution 
images  of  radar  targets  from  delay-doppler,  spotlight-mode  radar 
data.  This  approach  is  based  on  a  model  for  the  target's  reflectivity 
in  terms  of  wide-sense  stationary,  uncorrelated  scatterers  haying 
complex-valued  Gaussian  statistics.  The  imaging  problem  is  to 
estimate  the  target’s  scattering  function  in  terms  of  radar-echo 
data  acquired  with  a  series  of  target  illuminations.  We  develop 
a  method  for  solving  this  multidimensional  spectrum  estimation 
problem  through  the  use  of  maximum-likelihood  estimation 
implemented  via  the  expectation-maximization  algorithm. 


INTRODUCTION 

A  system  for  forming  the  image  of  a  radar  target  is  shown  in 
Fig.  I.  There  are  two  modes  for  collecting  data  to  form  the  image. 
An  antenna  of  sufficient  size  may  be  used  to  focus  radar  energy 


onto  a  patch  of  the  target  having  a  size  corresponding  to  a  resolution 
element  By  varying  the  position  of  the  incident  energy  in  some 
form  of  raster  pattern  over  the  target,  all  reflecting  patches  may 
be  illuminated  with  a  series  of  radar  pulses  and  an  image  of  the 
target  formed  by  displaying  the  return  energy  for  each  patch. 
Alternatively,  in  spotlight  mode,  the  energy  is  relatively  unfo¬ 
cused,  and  the  entire  target  illuminated  simultaneously  by  each 
transmitted  pulse.  The  same  rrm  ol  image  of  energy  versus  range 
and  cross-range  coordinates  can  be  formed  from  the  more 
complicated  echo  data  by  suitable  processing  which  utilizes  delay 
and  doppler-shift  sanations  present  in  a  scries  of  target  illumi¬ 
nations  Our  concern  is  wnh  forming  high  resolution  images  from 
data  acquired  in  the  spotlight  mode. 

There  are  at  least  two  uses  for  high  resolution  images  of  a 
radar  target.  One  is  in  developing  a  catalog  of  radar  cross-section 
profiles  for  various  target  types.  Another  is  for  target  identifi¬ 
cation.  The  latter  use.  illustrated  in  Fig.  2,  normally  proceeds  in 
two  separate  and  independent  steps.  First,  the  target’s  image  is 
formed,  and  then  features  of  the  target,  such  as  edges  and  textures, 
are  extracted  and  used  with  any  collateral  information  that  may 
be  available  to  identify  the  target.  As  illustrated  in  Fig.  3,  we 
are  also  interested  in  developing  a  more  coordinated  approach  to 
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Figure  2.  The  use  of  an  image  for  identification. 


Figure  3.  Coordinated  target  imaging  and  feature  extraction. 


forming  an  image  of  the  target  and  extracting  features  so  that 
these  two  processes  may  interact  constructively  so  as  to  improve 
both.  As  a  result  of  other  developments  in  our  laboratory  {4,5], 
we  believe  that  the  method  we  have  developed  for  forming  target 
images  is  ideally  suited  for  making  this  potentially  important 
extension. 

Wehner  [  1  ]  and  Mensa  [2]  describe  methods  for  forming  images 
of  radar  targets  from  spotlight  mode  data.  A  series  of  identical 
radar  pulses  (or  pulse  groups)  having  a  linear  FM  chirp  modulation 
is  used  to  illuminate  the  target.  The  echoes  are  processed  with  a 
two-dimensional  Fourier  transform  to  form  the  image.  This 
approach  is  based  on  an  intuitive,  deterministic  analysis  which 
results  in  accurate  target  images  under  high  signal-to-noise  ratio 
conditions.  Our  approach  differs  by  incorporating  a  statistical 
model  for  the  target’s  reflectivity,  accommodating  receiver  noise, 
and  in  using  a  statistical  estimation  approach  for  developing  a 
method  for  forming  the  image.  A  full  development  of  our  method 
is  contained  in  [3]. 

TARGET  MODEL 

We  model  the  complex  envelope  of  the  echo  data  according  to: 

r(t)-  |  s  T  ( f  —  x)y(f  -  x  .  T)dT  *  U'(<).  (I) 

where  iu(i)  is  a  white  Gaussian  noise  with  spectral  density  ,V„, 
s-tn  is  the  complex  envelope  of  the  transmitted  signal,  and  y(t.  t) 
is  the  reflectivity  at  time  t  of  all  reflecting  patches  at  two-way 
delay  r  ,  y(t.  r)  may  be  expressed  in  terms  of  all  reflectivities  at 


delay  x  according  to 

y(f.t)- Jc(/.T)e'2,,,d/.  (2) 

where  c(/.t)  is  the  reflectivity  of  all  the  points  on  the  target  at 
two-way  delay  x  which  introduce  a  Doppler  shift./  .  Since  targets 
are  of  finite  extent,  both  y(< .  t)  and  c(/ .  t)  are  zero  for  x  outside 
some  fixed  interval,  and  c(/.t)  is  zero  for  /  outside  a  fixed 
interval  because  the  Doppler  variable  is  equivalent  to  the  cross¬ 
range  coordinate  of  a  rotating  target. 

As  we  develop  in  detail  in  [3],  (I)  and  (2)  may  be  discretized 
into  the  matrix-vector  form 

r-f'c-w.  (3) 

where  the  superscript  "h”  denotes  the  Hermitian  transpose  oper¬ 
ation,  r it  1,  c II 1,  and  \v it  1  are  vectors  of  samples  of  r(( ),  c(/.  t) , 
and  iv(0,  respectively,  and  r  is  a  matrix  each  element  of  which 
is  the  product  of  a  sample  of  sT(t)  and  a  complex  exponential. 

We  have  adopted  a  diffuse-target  model  for  the  reflectivity. 
For  this,  the  target  is  assumed  to  be  comprised  of  uncorrelated 
scatterers  each  of  which  introduces  a  complex-valued,  zero  mean 
Gaussian  random  variable  as  a  multiplicative  factor  on  the  incident 
signal  reflected  by  it.  The  superposition  of  these  according  to  (1) 
and  (2)  results  in  the  radar  echo  data.  This  assumption  implies 
that  the  reflectivity  vector  chas  a  Gaussian  distribution  with  zero 
mean  and  'diagonal  covariance  matrix  I»F(cc*).  This 
covariance  consists  of  samples  of  the  scattering  function  of  the 
target,  which  is  the  power-spectrum  of  the  reflectivity  process 
y(r.-c).  For  our  approach,  the  reflectivity  is  a  two-dimensional 
Gaussian  process,  and  the  scattering  function  is  its  spectral 
intensity.  Further,  the  received  vector  r  has  a  Gaussian  distribution 
with  zero  mean  and  covariance  matrix 

k,  -  r'zr  *n„i.  (-») 

The  loglikelihood  of  the  data  r  may  be  expressed  in  terms  of  this 
covariance  matrix  according  to 

£(1<  ,  :r)  -  -  ln(def  K  , )  -  r*K  j1  r.  (5) 

IMAGING  PROBLEM 

Two  different  images  of  the  target  may  be  formed,  one  being 
an  estimate  of  the  target’s  reflectivity  c  and  the  other  being  an 
estimate  of  its  scattering  function  I,  both  images  being  displayed 
in  range  and  cross-range  coordinates.  A  unique  aspect  of  our 
method  is  that  it  produces  a  maximum-likelihood  estimate  of  Z 
and,  also,  a  conditional  mean  estimate  of  c  so  that  both  of  the 
possible  images  of  the  target  can  be  displayed  if  desired. 

Maximizing  the  loglikelihood  (5)  subject  to  the  constraint  that 
K,  must  be  of  the  parameterized  form  in  (4)  leads  to  a  trace 
condition  first  discussed  by  Burg,  Luenberger,  and  Wenger  [6] 
and,  in  the  present  context,  by  Snyder,  O'Sullivan,  and  Miller  [3]. 
While  this  trace  condition  in  principle  specifies  the  maximum- 
likelihood  estimate  of  the  scattering  function,  it  is  a  highly 
nonlinear  equation  with  no  closed  form  solution.  For  this  reason, 
we  have  in  [3]  adopted  the  use  of  the  alternating  maximization 
approach  of  Dempster,  Laird  and  Rubin  [7]  for  producing  the 
maximum-likelihood  estimate  numerically. 

A  sequence  of  estimates  of  the  scattering  function  is  obtained. 
The  estimate  at  step  p  *  1  is  defined  by  the  conditional  expectation 

I"’"’  -  ffcc*  1 11'”,  r  ].  (6) 

d 

where  ”  -  "  means  that  the  off  diagonal  terms  on  the  left  are  zero 
and  the  diagonal  terms  equal  the  diagonal  terms  on  the  right. 
Evaluating  the  right  side  of  (6)  is  a  standard  problem  in  estimation 
theory;  it  is  given  by 

i'"  ’♦  i(",r(r',z'',,r*  ,v0i)'' 

x(rr,‘-r"5:<'’T-Ivoi)(r',z<",r*,v0i)''r',i<'”.  (7) 

Equations  (6)  and  (7)  define  an  iteration  sequence  which  produces 
an  estimate  of  the  scattering  function  at  each  step.  Thus,  at  each 
step,  an  image  of  the  target  may  be  displayed,  and  each  image 
has  successively  higher  likelihood.  At  each  step,  the  conditional 
mean  estimate  of  the  reflectance  is  also  generated.  At  step  p.  this 
estimate  is 

ctriz''”.  rj- i<'”r(r,z,'"r  *  ,v„i)  v  (s> 

EXAMPLE 

At  the  present  time,  we  are  implementing  (6)  and  (7)  to  run 
on  a  Distributed  Array  Processor  (DAP)  made  by  the  Active 
Memory  Technology  Company,  which  is  a  mesh  connected  parallel 


computer  with  1024  processors.  This  implementation  will  permit 
comparisons  to  be  made  between  this  maximum-likelihood  method 
and  the  more  conventional  method  that  employs  two-dimensional 
Fourier  transforms. 

We  have  performed  a  preliminary  computer-simulation  study 
in  which  there  is  a  point  target  that  is  concentrated  at  a  single 
range  and  crossrange.  For  this  situation,  forming  the  target  image 
corresponds  to  estimation  of  the  power  spectrum  of  a  time  series 
having  one  nonzero  spectral  component.  Fig.  4  shows  the  result. 
The  graph  on  the  left  in  Fig.  4  is  the  output 
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Figure  4.  Shown  are  the  output  SNR  (left)  and  bias  (right) 
for  estimating  the  reflected  power  of  a  point  target  in  noise. 


signal-to-noise  ratio  versus  the  input  signal-to-noise  ratio,  and 
the  graph  on  the  right  is  bias  versus  the  input  signal-to-noise 
ratio,  where  these  quantities  were  estimated  from  3000  inde¬ 
pendent  trials  and  are  defined  according  to 

SNRln-F0/A/0,  (9) 

where  Zf0  is  the  average  power  in  each  of  the  P  nonzero  spectral 
components  of  the  reflectivity,  as  defined  by 

Fo-  jjTr(I): 


and 


SNR. 


Z 

v  MSE(£) 


(10) 


where  MSEC  £ )  is  the  sample  mean-square  error  in  estimating  E; 
and 

BIAS-S-I.  (ll) 


The  estimates  of  the  scattering  function  for  the  maximum- 
,  likelihood  method  and  for  the  two-dimensional  Fourier  transform 
method  (which  in  one  dimension  becomes  a  periodogram)  are 
compared  in  Fig.  4.  The  maximum-likelihood  estimates  (new 
method)  are  seen  to  have  high  output  SNR  and  low  bias  compared 
to  the  periodogram  for  input  SNRs  above  5  dB.  This  very  pre¬ 
liminary  result  gives  us  optimism  that  superior  images  of  diffuse 
targets  will  be  obtained  with  the  maximum-likelihood  method. 

RELATED  ISSUES 

In  [8],  we  have  developed  a  Cramer-Rao  lower  bound  on  the 
mean-square  error  performance  of  maximum-likelihood  estimates 
of  the  scattering  function.  The  Fisher  information  matrix  of  this 
bound  is  shown  in  [8]  to  play  an  important  role  in  specifying 
conditions  for  the  uniqueness  of  estimates  of  the  scattering 
function  and  in  the  selection  of  radar  pulse  shapes  for  achieving 
low  variance  images.  Convergence  of  the  iterative  sequence  in 
(6)  and  (7)  is  discussed  in  (9).  In  [10],  we  develop  a  method  for 
incorporating  equality  and  inequality  constraints  on  the  input 
signal-to-noise  ratio,  which  results  in  improved  estimates  when 
such  side  information  is  available. 
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ABSTRACT 

Maximum-likelihood  spectrum  estimation  is  an  ill-posed  problem.  In  this  paper,  we  use  a  method  of 
sieves  for  addressing  this  issue.  The  estimate  of  the  spectrum  is  constrained  to  a  subset  of  some  Hilbert 
space  of  functions  over  which  a  complete  set  of  nonorthogonal  basis  functions  is  defined.  The  estimate  is 
then  represented  by  a  countable  set  of  coefficients  in  a  nonorthogonal  series  expansion.  By  defining  an 
appropriate  sieve  on  this  countable  set,  our  problem  reduces  to  maximum-likelihood  estimation  of  the 
parameters  in  the  sieve.  Three  main  attractive  features  of  this  approach  are:  (1)  the  nonorthogonal  expan¬ 
sion  is  a  convenient  framework  for  defining  the  sieve  and  including  a  priori  information;  (2)  mean-square 
consistency  of  the  estimates  can  be  expected;  and  (3)  we  have  derived  a  tractable  alternating  maximization 
algorithm  for  estimating  the  parameters.  The  setup  of  this  problem  is  general  and  can  be  applied  without 
major  difficulties  to  the  estimation  of  higher-dimensional  spectral  functions,  as  occurs,  for  example,  in 
imaging  radar  targets  from  delay-doppler  data. 

1.  INTRODUCTION 

Maximum-likelihood  (ML)  estimation  of  a  continuous  function,  or  even  of  an  infinitely  countable  set 
of  parameters  is  known  to  be  an  ill-posed  problem  [1],  In  this  paper,  we  use  a  regularization  method  for 
addressing  this  issue  when  the  continuous  function  is  the  spectral  density  of  a  Gaussian  process,  and  this 
density  must  be  estimated  from  one  single  sample  function  of  the  process.  The  ill-posedness  of  the  estima¬ 
tion  problem  is  often  addressed  by  approximating  the  model  for  the  process  by  a  simpler  one,  assuming 
temporal  periodicity.  The  spectral  density  is  then  discrete.  It  follows  that  for  bandlimited  processes,  the. 
spectrum  is  finite-dimensional  and  can  be  estimated  by  means  of  standard  ML  techniques.  However,  the 
fundamental  difficulties  inherent  to  ill-posed  problems  remain,  indicating  that  the  assumption  of  a  periodic 
process  is  not  sufficient.  For  instance,  the  length  of  the  periodic  extension  of  the  data  cannot  be  made  arbi¬ 
trarily  large  without  the  estimates  becoming  very  rough.  Furthermore,  the  estimate  obtained  from  the 
observation  of  one  single  realization  of  the  process  is  not  consistent,  no  matter  how  many  samples  are  col¬ 
lected. 

The  approach  we  use  for  estimating  the  continuous  spectral  density  does  not  require  any  approxima¬ 
tion  of  the  model.  It  is  based  upon  Grenander’s  method  of  sieves,  which  provides  a  framework  for  estima¬ 
tion  in  an  infinite-dimensional  parameter  space  [1].  The  estimation  is  performed  in  a  restricted  parameter 
set  over  which  the  estimates  are  stable.  For  an  infinite  amount  of  data,  this  restricted  parameter  set  is  dense 
in  the  actual  parameter  set.  This  procedure  leads  to  producing  stable  and,  we  expect,  consistent  estimates. 

Section  2  of  this  paper  contains  a  description  of  the  model  for  the  stochastic  process  whose  spectral 
density  is  sought.  In  Section  3,  we  briefly  summarize  some  previous  results  in  the  literature,  in  which  an 
approximation  of  the  model  is  made;  the  limitations  of  these  methods  are  mentioned.  Section  4  contains  a 
short  review  of  regularization  methods.  The  sieve-constrained  ML  estimator  is  described  in  Section  5. 

•  The  work  described  in  this  paper  was  supported  by  contract  number  N00014-86-K-0370  from  the  Office  of  Navtl 
Research  and  by  grant  number  MIP-8722463  from  the  iVaaonai  Science  Foundation. 
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2.  MATHEMATICAL  MODEL 

Consider  a  complex-valued,  wide-sense  stationary  stochastic-process  y(t).  Our  estimation  pro¬ 
cedure  is  based  upon  the  following  statistical  model  for  the  observations.  The  process  is  Gaussian  with  a 
mean  of  zero.  Its  spectrum  S  (f )  is  bandlimited  to  some/maI,  and  N  samples  of  one  realization  of  the  pro¬ 
cess  are  observed,  with  the  sampling  frequency  1/A/  at  least  equal  to  the  Nyquist  rate  2 Furthermore, 
the  process  is  observed  in  complex-valued  additive  white  Gaussian  noise  w(r )  with  spectral  density  A'q. 

Under  the  assumptions  stated  above,  the  samples  of  y  (/)  can  be  represented  via  the  Cramer  spectral 
representation: 

/- 

y(nA/)  =  j  exp(-j2KfnAt)Z(df)  ,0<n<N-l.  (2.1) 

In  this  generalized  Stieltjes  integral,  Z(df)  is  an  orthogonal,  Gaussian  spectral  process  with  variance 
S(f)df[  l,p.74]. 

In  the  next  sections,  a  ML  estimator  for  S(f )  is  presented.  We  shall  find  it  convenient  to  consider  a 
discrete  approximation  to  the  integral  (2.1).  Define  a  uniform  partitioning  of  the  interval  [-/^,,/m,,], 

(/* }  t  ,  with  partition  size  Aw  =  2/^JM .  Next  define  the  process 

c(k)  =  (Aw)~l/2  [Z(f i„.i)-Z(f  *)]  ,0<k<M-\.  (2.2) 

This  process  is  Gaussian  and  orthogonal.  The  variance  of  each  sample, 

/*•> 

02(k)  =  -7—  (2.3) 

is  an  approximation  to  S(fk).  Next,  approximating  the  Stieltjes  integral  (2.1)  with  a  Riemann  sum,  we  get 

/(«AO«  V -^T^^'exp  (-jlnfknAt)c  {k) .  (2.4) 

M  i=0 

As  the  partition  size  tends  to  zero,  y  (.)  converges  to  the  stochastic  integral  y  (.)  in  the  mean-square  sense: 

y(nAt)  =  l.i.m.y(nAt)  (2.5) 

M  — *  •» 

Notice  that  in  the  special  case  of  a  periodic  process  with  period  M  At,  y  (.)  =  y  (.),  so  that  (2.4)  is  the 
exact  spectral  representation  of  y  (n  At).  This  observation  has  motivated  using  (2.4)  as  a  model  for  y  (.)  in 
[2].  The  parameter  M  is  greater  or  equal  to  N  but  is  not  further  specified,  but  Af  should  be  large  enough  to 
ensure  that  the  model  approximation  is  valid.  In  [2]  the  issue  of  ill-posedness  is  addressed  by  keeping  M 
finite.  In  the  next  section,  we  will  summarize  some  known  results  for  this  estimator.  On  the  other  hand, 
the  sieve-constrained  ML  estimator  presented  here  is  still  based  on  the  representation  (2.4),  but  does  not 
require  Af  to  be  finite.  As  such,  this  estimator  does  not  involve  any  model  approximation,  as  shown  in 
(2-5). 

We  now  specify  the  model  for  the  observations.  Using  the  same  vector  notations  as  in  [2],  we  get, 
from  (2.4): 

r  =  +  w  ,  (2.6) 

where  r  is  the  N- vector  of  observations,  w  is  a  noise  vector  with  diagonal  covariance  N  ,  where  IN  is  the 
NxN  identity  matrix,  and  c  is  the  M-vector  formed  from  the  samples  of  the  orthogonal  process  c(k) 

defined  in  (2.2).  Since  c  is  orthogonal,  its  covariance  I  is  diagonal  with  entries  (^(k))  .  T  is  a  MxN 

matrix  with  (k,n)  entry  given  by  \  -~^-exp(-j2KfknAt):  finally,  t  denotes  the  Hermitian-transpose 
operator  for  matrices  * . 


T  can  be  generalized  readily  to  accomodate  a  linear  tranjformauon  of  y  (.)  in  the  measurement,  as  occurs  in  the  radar 

-  43  - 


-3- 


The  covariance  matrix  for  r  is  given  by 

Kr  =E[rrf]  =  I^EF  +  Nq!s  .  (2.7) 

3.  UNCONSTRAINED  ML  ESTIMATOR 

In  this  section,  we  introduce  the  ML  estimator  for  the  spectral  matrix  2  defined  in  the  model  (2.6). 
From  (2.7),  the  loglikelihood  function  for  2  is 

L(r  ,Z)  =  -  \n  deUX'lT  +  NqI")  - rWzr  +  NoIsTlr  .  (3.1) 

Maximizing  the  likelihood  with  respect  to  2  yields  the  necessary  trace  condition  which  a  positive  definite 
estimate  2  must  satisfy  [2]: 

Tr  [TA7V -^)A7!r'S2]  =  0  ,  (3.2) 

with 

for  all  admissible  M  xAf  diagonal  matrices  82  [2].  This  trace  condition  is  a  nonlinear  equation  in  2.  Gen¬ 
erally  it  cannot  be  solved  directly  in  closed-form,  so  some  numerical  search  procedure  must  be  imple¬ 
mented.  An  approach  is  the  expectation-maximization  (EM)  algorithm  of  Dempster  et  al.  [3]  used  in  [2]. 
The  complete  data  are  defined  as  (c,w).  An  initial  estimate  2(0)  is  selected.  At  step  k+1  (k  =  0,1,..)  the 
estimate  is  updated  according  to 

=  argmax  Q  (2|2(i))  (3.3) 

where 

Q  (2|2(*>)  =  -  In  o*(i )  -  "£  ,  (3.4) 

;=o  i=o  cr(i) 

and 

E[|c(i)|2|r,2<t)]  =  [i(k)  -  0-k)rKr-1  ^rri(k)  +  i{k)rK~l (k)rr fK~l (i)rf2(i)]  (jj) ,  (3.5) 

where 

Krik) = . 

This  algorithm  produces  a  sequence  of  estimates 

62(i)(*+1)  =  E[|c(i)|2|r,Z<*)l  (3.6) 

having  increasing  likelihood.  It  can  be  shown  that  the  stable  points  of  this  algorithm  satisfy  the  necessary 
trace  condition  (3.2)  for  a  maximizer  [2], 

The  performance  of  this  estimator  has  been  evaluated  in  previous  studies  and  compared  to  traditional 
spectrum  estimators  such  as  the  periodogram  [4, 5],  It  has  been  shown  that 

(1)  the  ML  estimates  have  a  very  low  bias; 

(2)  the  ML  estimates  offer  a  better  trade-off  between  variance  and  resolution  than  the  periodogram; 

(3)  the  ML  estimates  are  not  mean-square  consistent  This  unfortunate  result  is  due  to  the  fact  that  the 
dimension  M  of  the  parameter  space  is  at  least  equal  to  the  number  of  data  N .  We  are  actually  deal¬ 
ing  with  a  small-sample  ML  estimation  problem.  This  observation  is  one  factor  that  motivated  the 
study  presented  in  the  next  sections. 


imaging  problem  {2). 
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4.  REGULARIZATION  OF  THE  ML  ESTIMATES 


According  to  Hadamard’s  classical  definition,  a  problem  is  said  to  be  well  posed  if:  (1)  it  has  a  solu¬ 
tion,  (2)  the  solution  is  unique,  and  (3)  the  solution  varies  continuously  with  the  data.  A  problem  is  ill 
posed  if  it  is  not  well  posed.  Small  changes  in  the  data  can  then  produce  unbounded  changes  in  the  esti¬ 
mates.  Typical  examples  are  inversion  of  certain  Fredholm  integral  equations  and  the  reconstruction  of 
functions  from  truncated  Fourier  transforms  [6]. 

The  ML  estimation  procedure  produces  similar  artifacts  when  the  number  of  parameters  (M)  is  large 
compared  to  the  number  of  data  samples  ( N ).  The  estimates  are  extremely  sensitive  to  small  changes  in 
the  data  and  exhibit  a  very  rough  shape  characterized  by  very  sharp  peaks  and  low  valleys.  We  can  say 
that  the  problem  is  practically,  if  not  technically,  ill  posed.  Such  a  behavior  is  highly  undesirable,  and  a 
regularization  of  the  estimates  is  required.  This  concept  was  introduced  by  Tikhonov  [7].  A  large  bulk  of 
mathematical  literature  has  been  written  on  this  subject;  see  for  instance  Bertero  [8]  for  a  tutorial. 

One  possible  method  for  the  regularization  of  estimates  is  the  method  of  sieves,  introduced  by 
Grenander  [1].  According  to  Grenander’s  definition,  a  sieve  in  a  parameter  space  A  is  a  family  of  subsets 
S(y)  of  A  indexed  by  a  positive  parameter  y  called  the  mesh  size.  A  restricted  ML  estimate  exists  over 
each  set  S(y).  As  the  mesh  size  tends  to  zero  the  sets  S(y)  will  be  large  enough  to  allow  the  ML  solution  to 
converge  to  any  solution  in  A  [1,  p.357].  For  the  problem  at  hand,  the  choice  of  a  "good"  value  of  y  will 
depend  on  the  data  record  size,  the  noise  level,  and  possibly  other  factors.  A  major  problem  is  to  find 
sieves  which  will  make  the  ML  estimates  converge  and  possess  some  desirable  practical  features  such  as 
analytical  and/or  computational  tractability.  For  instance,  one  might  consider  a  convolution  sieve: 

S(y)=j  cr(t)  |  cr(i)  =  T|(/)  *  VjxO)  o <Af-l| ,  (4-1) 

where  cr(i)  is  the  constrained  spectrum,  T](i)  is  any  spectrum,  and  ^(t )  is  a  convolution  kernel,  or  "win¬ 
dow”,  with  mesh  equal  to  y.  However,  it  is  not  clear  that  finding  spectral  estimates  within  this  sieve  is 
computationally  tractable.  In  1985  Chow  and  Grenander  recognized  that  regularizing  the  ML  spectrum 
estimates  with  a  convolution  sieve  appears  to  be  a  formidable  task,  so  they  studied  other  types  of  estima¬ 
tors  [9]. 

Here  we  recommend  a  regularization  method  based  on  nonorthogonal  expansions  of  the  unknown 
density  function.  This  approach  leads  to  computationally  tractable  ML  estimates.  As  discussed  in  §5.3,  we 
also  expect  the  estimates  to  be  asymptotically  consistent 

The  regularization  procedure  developed  in  the  next  section  will  be  applied  to  the  technically  ill-posed 
spectrum  estimation  problem  in  which  the  ratio  M  IN  tends  to  infinity,  where  the  approximate  model  (2.4) 
converged  to  the  exact  model  (2.1).  We  will  show  how  to  solve  this  infinite-dimensional  parameter  estima¬ 
tion  problem. 

5.  A  METHOD  OF  SIEVES  BASED  ON  NONORTHOGONAL  EXPANSIONS 

The  spectrum  S  If)  belongs  to  a  Hilbert  space  of  functions  defined  over  [-/„„,/  J.  Denote  by  H 
this  Hilbert  space,  and  by  H,.  the  subset  of  all  nonnegative  functions  in  H.  Typically  H  is 
L2  space  of  square-integrable  functions  over  Our  goal  is  to  represent  the 

spectral  function  in  H+  in  terms  of  a  set  of  basis  functions  having  localized  support.  This  permits  a  good 
local  representation  of  a  function.  As  such,  this  is  related  to  the  concept  of  wavelets. 

5.1.  Definition  of  the  sieve 

Consider  a  spectral  M -vector  or.  To  this  vector  we  associate  a  step  function  S(f)  according  to  the 
map  S7W  defined  as  follows; 

ST^to2]  :  R *  H  :  <T  -»  S(/)  =  ST„  [a2]  ,  (5.1) 
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M-l 

where  {/,■ }  is  the  same  uniform  partition  of  the  interval  [-/m„./ra„]  as  in  Section  2,  RM  is  the  M- 
dimensional  Euclidean  space,  and  Xiif)  *s  the  indicator  function  over  the  i-th  partition  interval. 

Let  (Vm(/))m=o  be  a  nonorthogonal  set  of  positive  basis  functions  in  H+.  Consider  the  set  A  of  all 
functions  in  the  span  of  the  (ym  ]  with  positive  coefficients: 

A=js(/) \S(f)=  £a(m)Vm<y)  ,a(m)>oJ.  (5.2) 

We  define  our  parameter  set  to  be  A.  Clearly  A  is  a  subset  of  H+.  It  would  be  desirable  to  have  (ym  ] 

be  designed  so  that  A  is  dense  in  H+.  It  is  well  known  that  this  can  be  achieved  with  defined  as 

indicator  functions  over  a  partition  of  the  frequency  domain,  in  the  limit  as  the  partition  size  tends  to  zero. 
How  A  is  made  dense  in  H+  for  more  elaborate  designs  remains  to  be  investigated. 

Now  define  the  following  step-function  approximation  of  \| r„(f): 

(5.3) 

where  \j/«^  is  a  M  -vector  of  samples  of  ym(/)  taken  on  the  partition  of  the  frequency  interval.  As 
M  ¥**,*,(/)  converges  to  (f )  in  the  norm  of  H.  We  define  a  sieve  S(M ,Q)  on  S(f)  in  A  by  trun¬ 
cating  the  series  representation  (5.2)  to  a  finite  number  of  terms  Q  and  considering  M-step  approximations 
to  the  basis  functions, 

S(Af, Q)=|  S(f) \S(f)*  Q£a{m)yM„(f)  ,a(«)>o[.  (5.4) 

t  m=0  J 

The  functions  in  the  set  S(A { ,Q)  are  constrained  to  be  M-step  functions  of  the  type  (5.1).  For  the  time 

being,  we  keep  M  finite,  but  ultimately  M  will  be  allowed  to  tend  to  infinity.  We  view  as 

being  Q  unknown  parameters  for  which  ML  estimates  are  sought. 

Why  nonorthogonal  basis  functions? 

Before  considering, the  general  problem,  we  mention  the  special  case  where  the  basis  functions  con¬ 
sidered  in  (5.4)  are  Q  indicator  functions  over  disjoint  intervals  covering  the  frequency  domain.  In  this 
instance,  the  computation  of  the  a(m)’s  is  straightforward  *.  However  these  basis  functions  have  the 
major  disadvantage  that  the  estimates  are  step  functions,  and  so  display  discontinuities  that  are  uncharac¬ 
teristic  of  spectral  densities. 

To  achieve  good  representations  for  a  richer  class  of  functions,  we  have  to  forgo  the  aforementioned 
restriction  on  the  structure  of  the  basis  functions.  A  possible  design  of  the  basis  functions  might  be  based 
upon  dilations  and  translations  of  a  smooth  basic  function.  Such  a  design  is  used  in  the  wavelet  literature 
because  it  provides  a  convenient  multiscale  representation  of  the  function  of  interest  [10,  11].  This  is  usu¬ 
ally  achieved  more  easily  with  nonorthogonal  basis  functions  [11]. 

Our  strategy  in  using  other  basis  functions  will  have  practical  interest  if  a  tractable  parameter  estima¬ 
tion  procedure  can  be  derived.  This  issue  is  addressed  in  the  next  section. 

52  Computation  of  the  sieve  estimator 

To  each  function  S(f)  in  the  constrained  set  (5.4)  is  associated  a  spectral  M -vector  o2  via  the 
inverse  map 

=  STf  lS(f)] .  (5.5) 

•  The  maximization  problem  (5.7)  can  then  be  solved  in  dosed  form. 
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Since  the  map  ST^1  is  linear,  we  conclude  from  (5.1),  (5.4)  and  (5.5)  that  the  vector  ct2  is  constrained  as 
follows. 

,  e-t 

£2=  2>(m)\£M,,„  •  (5.6) 

m=0 


The  ML  estimates  of  the  a(m)'s  can  be  obtained  by  using  the  EM  algorithm  described  in  Section  3.  At 
each  step  the  expression  (3.4)  is  maximized  with  respect  to  ct2  subject  to  the  sieve  constraint  (5.5).  This 
amounts  to  finding  the  maximum  with  respect  to  a  of 


Q(a\d(k))  =  -  "£  In 

i=0 


a- 1 

£  a  («)>£»,*(/) 

m=& 


E[\c{i)\2\r4W ] 
i=0  2fl(*)v^(/) 

m=0 


(5.7) 


This  maximization  problem  is  very  difficult  and  must  be  solved  via  some  numerical  procedure.  Keeping  in 
mind  that  this  operation  has  to  be  repeated  at  every  iteration  of  the  EM-algorithm,  one  is  hardly  attracted 
by  this  prospect.  We  propose  the  following  method  instead. 

.An  EM  algorithm  has  been  developed,  based  on  the  definition  of  a  new  complete  data  space  for  pro¬ 
ducing  the  estimates.  This  method  is  set  up  as  an  unconstrained  maximization  problem  in  a  Q  -dimensional 
parameter  space  and,  as  such,  it  does  not  involve  any  computation  in  H.  Because  of  the  existence  of  such 
an  algorithm,  our  estimation  procedure  based  on  nonorthogonal  expansions  turns  out  to  be  computationally 
feasible.  This  usually  represents  an  impressive  computational  saving  and  is  a  major  attractive  feature  of 
this  approach. 


Complete/Incomplete  Data  Spaces 


Write  c  as  the  sum  of  Q  independent  vectors: 

c  =Q£cm,  (5.8) 

/n=0 

where  cm  is  a  M-vector,  sample  of  a  0-mean  Gaussian  process  with  diagonal  covariance  a  (m  )vi'(n ,  with 
a  diagonal  matrix  made  of  {  yMi*,(t),  i  e  [0..A/-1]  }.  These  Q  processes  are  independent.  Now 
define  the  Jm -vector  c„  as  the  restriction  of  cm  to  its  support  set  S„ .  The  covariance  '¥m  of  cm  is  made  of 
the  nonzero  entries  of  'f'm.  With  these  notations  our  model  (2.6)  becomes 
(2-i  . 

r  =  2  r  C‘m  +  w 
«=0 

=Q£rLcm+w .  (5.9) 

m=C 


We  define  the  complete  data  as  (  (c„  jSi  ;  w).  Following  this  definition  we  write  the  loglikelihood 
for  the  complete  data  as 
(2-1 

(<*(*)  =  Llnp(cm:a) 

mm 0 

=  -  Z  In  dot  (fl(«)VJ- cl  (a  (m  )Vm  )~l  cB  . 

mm 0 


Discarding  all  terms  not  involving  a ,  we  obtain: 


MOI2 


£2-1  S-i  1 

lcAa)  =  ~  ZJ*to*0n)-  ITT" T  L  - - 7J7  ■ 

The  conditional  expectation  of  l;d(a)  is  then 

Q  ( a  \d(i))  =  - 2£  JJna{m)~  V  -j—  V 

m-Oa(mhsS. 


E[lc,(i)|2|r.a(t)] 


(5.10) 


(5.11) 
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Taking  the  partial  derivative  of  (5.1 1)  with  respect  to  a  (m ),  and  setting  the  result  to  zero  yields 

--(»>)«•"- j-£  .  (mi 

3 n  ies.  VMjhO) 

Next  we  evaluate  the  conditional  expectation  of|cm(i)|2  in  (5.12).  Define 

c™  =  E[cn\r4W}.  (5.13) 

The  conditional  expectation  in  (5.12)  is  the  i-th  diagonal  element  of  the  matrix 
E[cncl\r4(k)]  =  £[ (cn-cik))  {cn-cik)Y\r,d(k)}  + 

=  (K™ Krl  (k)  r WY} KY  -  K^)Kr’lwK^  ,  (5.14) 

where  we  denote  by  Kv  the  conditional  correlation  E  [xyf\  d(i)]  of  two  random  vectors  x  and  y .  Now,  the 
expectations  are  evaluated  from  the  model  (5.9).  Note  that  because  the  cm ’s  are  independent,  the  condi¬ 
tional  expectations  on  d{k)  are  simply  conditional  expectations  on  a(m)(k\  We  get 

K?J  =  £  [  cy\d(k)]  =  d(m  rm  ,  (5.15a) 

=  E  ( rr'|d(i)]  =  £,a'0')(*)r/'P;  T,  +  Nq[n  ,  (5.15b) 

i- o 

=E  [cmcj,\dw]  =  d(mf)'¥„  .  (5.15c) 

Substituting  the  expressions  (5.15a,c)  into  (5.14)  yields 

E[cncl\r  ,d(k)]  =  d (m )w'Vn  +  d ( m  Tn  K'1  (k)[rr f-K^ )  K~x  (i)  d (m )(i) ,  (5. 16) 

with  Kr(k)  given  in  (5.15b).  Taking  the  i-th  diagonal  element  of  (5.16)  and  substituting  in  (5.12),  we  get 

d(m )(*+l)  =  7-  I  --1  a(m)(*VAf/n(0  +  d("O(tV«,«0‘) 

Jn  i«S„  VMjnV)  L  —  ~ 

x[r„  K~'  (k\rrf-Kr(k))Kr-Hk)ri  ]„■  VM^,(Oa(«)(t)] 

=  d (m )(t)  +  d (m )(i)  y-£^(0 

•  J  m  i€Sm 

xtr,  ]g  ]  d(m)ik).  (5.17) 


Now  defining  the  NxN  matrices 


Z(k)  =  K~x  {k\rrf-KXky)K-1  (i) 


k  —  —  -  rfvi/  v 

lxm  t  *  m  T*» 1  /*  * 
^  m 


the  term  between  brackets  in  the  right-hand  side  of  (5.17)  can  be  written 
7-  I  v.^0')(rmz<‘)rj;),  =  j-Tr  [  ^rMz(i’r: 

=  rr  J-rJXr^ 

J  m 

=  7>[a’„z(*)  . 

Substituting  (5.20)  in  (5.17)  we  get  the  update  equations  at  stage  k  of  the  algorithm: 

d(mf+l)  =  d(m)(k)  +  d(mfk)Tr  K„Zik)  d(m)ik) , 


-  48  - 


-8- 


Kr^  =  Q-£d{jYk)Kj  +NaIs  . 
i=o 


(5.22) 


Comments 

(1)  The  update  equations  (5.21)  and  (5.22)  highlight  the  role  played  by  the  basis  covariances  Km  in  this 
representation. 

(2)  As  M  the  algorithm  for  estimating  the  coefficients  a(m)  is  well-behaved.  The  reason  is  that 
the  basis  covariance  Kn  in  (5.21)  is  then  given  by  the  limit  of  the  matrix  product  (5.19)  as  M  -» 
which  is  a  matrix  of  integrals  with  entry  (n  jn)  given  by 

/- 

Km{n,rn)=  j  exp (-/ 2nf (n-m)At)  \ym(f)df  . 

We  can  thus  let  M  «>  without  any  loss  in  stability  of  the  estimates  or  any  increse  in  computations. 

(3)  Implementation  of  the  solution: 

Each  matrix  Km  is  constant  and  can  be  computed  off-line.  At  each  step  the  covariance  matrix  Kr  in 
(5.22)  must  first  be  inverted,  which  accounts  for.V3  operations.  Following  this,  updating  each  aim) 
just  requires  N 2  multiplications/additions.  The  complexity  of  each  step  of  the  EM  algorithm  is  then 
N3  +  QN 2,  as  compared  to  /V3  +  MN1  in  the  unconstrained  case.  The  saving  is  impressive  for  infinite 
M,  keeping  Q  ~N. 

(4)  It  is  easily  shown  that  the  trace  appearing  in  the  update  equation  (5.21)  is  the  partial  derivative  of  the 
incomplete-data  loglikelihood  with  respect  to  a  ( m ). 

5  J.  Some  Regularization  Aspects 

The  most  attractive  feature  of  this  approach  is  a  theoretical  one.  We  are  effectively  estimating  Q 
parameters  instead  of  M  as  in  the  unconstrained  problem.  The  choice  of  Q  is  somewhat  arbitrary  and 
depends  on  the  way  we  have  decided  to  represent  the  spectral  function.  As  a  rule  of  thumb  it  seems  desir¬ 
able  to  have  Q  (#  of  parameters)  no  larger  than  N  (#  of  data).  This  can  be  done  while  letting  M  be  infinite, 
a  result  that  could  not  have  been  obtained  without  regularization  of  the  ML  estimates.  The  process  in  our 
model  then  approaches  the  limiting  case  of  a  (nonperiodic)  stationary  process,  and  as  such  does  not  involve 
any  model  approximation. 

Consistency  of  the  estimates 

A  major  goal  when  defining  a  sieve  is  to  ensure  consistency  of  the  estimates.  Although  no  thorough 
study  has  yet  been  undertaken  for  the  sieve  (5.4),  we  believe  that  consistency  can  be  achieved,  provided 
that  Q  tends  to  infinity  at  a  slower  rate  than  N.  Typically  :  Q  -  Na  ,  with  0  <  a  <  1.  Another  objective  is 
to  ensure  rapid  convergence  of  the  estimates.  Clearly  the  issue  here  is  the  design  of  the  basis  functions, 
which  is  one  of  our  current  research  areas. 


CONCLUSIONS 

We  have  proposed  a  method  for  a  estimating  a  spectral  function  represented  by  a  linear  combination 
of  basis  functions.  The  task  of  finding  the  ML  estimates  of  the  coefficients  of  these  basis  functions  looks 
formidable  at  first  (refer  to  the  optimization  problem  (5.7)).  However  we  propose  a  method  based  on  the 
decomposition  (5.8)  of  the  spectral  process  into  a  sum  of  Q  independent  spectral  processes  with  known 
covariance,  up  to  a  scale  factor.  The  subsequent  reformulation  of  the  complete/incomplete  data  spaces 
leads  to  the  derivation  of  a  powerful  algorithm  for  evaluating  these  scale  factors. 
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The  main  features  of  this  approach  are:  (1)  flexibility.  The  choice  of  basis  functions  is  wide  open.  In 
particular  we  can  choose  the  shape  of  the  basis  functions,  amount  of  overlap,  and  support  set  extent  (possi¬ 
bly  non-uniform).  (2)  using  an  exact  model.  No  approximation  of  the  exact  model  for  the  stochastic  pro¬ 
cess  is  needed.  (3)  computational  efficiency.  The  computational  complexity  of  the  algorithm  is  a  function 
of  the  dimension  of  the  parameter  space,  not  of  the  partition  size  of  the  frequency  axis.  An  important  result 
is  that  it  is  now  possible  to  let  M  — » =»  without  any  increase  in  the  number  of  compn'adons 

Several  interesting  questions  remain  open  issues.  In  particular  it  is  desircu  .j  ucmeve  mean-square 
consistency  of  the  estimates.  As  mentioned  in  §5.2.2,  this  might  be  done  by  letting  the  sieve  grow  at  an 
appropriate  rate  as  jV  — » The  issue  of  uniqueness  of  the  estimates  should  also  be  investigated.  Finally, 
the  design  of  the  basis  functions  will  determine  the  overall  estimator  performance,  as  far  as  convergence 
rate  and  sensitivity  to  noise  are  concerned. 

Also  note  that  the  setup  of  this  problem  is  very  general  and  can  be  applied  without  major  difficulties 
to  the  estimation  of  higher-dimensional  spectral  functions.  In  particular,  the  spectra  of  interest  in  the 
radar-imaging  problem  described  in  [2]  are  two-dimensional.  The  application  of  the  ideas  in  this  paper  to 
the  radar-imaging  problem  is  another  area  of  active  research. 

REFERENCES 

[1]  U.  Grenander,  Abstract  Inference,  Wiley  1981. 

[2]  D.  L.  Snyder,  J.  A.  O’Sullivan  and  M.  I.  Miller,  "The  Use  of  Maximum-Likelihood  Estimation  for 
Forming  Images  of  Diffuse  Radar-Targets  from  Delay-Doppler  Data",  Vol.  35,  No.  3,  IEEE  Trans, 
on  Information  Theory,  pp.  536-548,  1989. 

[3]  A.D.  Dempster,  N.M.  Laird  and  D.B.  Rubin,  "Maximum-Likelihood  from  Incomplete  Data  via  the 
EM  Algorithm",  J.  Royal  Scat.  Soc.,  Vol.  B39,  pp.  1-37,  1977. 

[4]  P.  Moulin,  "Maximum-Likelihood  Spectrum  Estimation  of  Periodic  Processes  from  Noisy  Data", 
Washington  University,  ESSRL  internal  report. 

[5]  P.  Moulin,  D.L.  Snyder  and  J.A.  O’Sullivan,  "Maximum-Likelihood  Spectrum  Estimation  of 
Periodic  Processes  from  Noisy  Data",  Proc.  of  the  1989  Conference  on  Information  and  Systems  Sci¬ 
ence,  Johns  Hopkins  University,  Baltimore,  1989. 

[6]  L.S.  Joyce  and  W.L.  Root,  "Precision  Bounds  in  Superresolution  Processing",  J.  Opt.  Soc.  Am., 
Void,  pp.  149-168,  1984. 

[7]  A.  Tikhonov  and  V.  Arsenin,  Solutions  of  Ill-Posed  Problems,  Winston  &  Sons,  Washington  D.C. 
1977. 

[8]  M.  Bertero,  C.  De  Mol,  G.  Viano,  "The  Stability  of  Inverse  Problems",  Chapter  5  in  Inverse  Scatter¬ 
ing  Problems,  Springer-Verlag  1980. 

[9]  Y.  Chow  and  U.  Grenander,  "A  Sieve  Method  for  the  Spectral  Density”,  Ann.  Star.,  Vol.  13,  No.  3, 
pp.  998-1010, 1985. 

[10]  M.  Frasier  and  B.  Jawerth,  "Decomposition  of  Besov  Spaces",  Indiana  Univ.  Math.  Journal,  Vol.  34, 
No.  3,  1985,  pp.777-799. 

[11]  I.  Daubechies,  A.  Grossmann  and  Y.  Meyer,  "Painless  Nonorthogonal  Expansions",  J.  Math.  Phys., 
Vol.  27,  No.  5,  1986. 


-  50  - 


6.5  Reprint  of:  J.  A.  O’Sullivan,  P.  Moulin,  D.  L.  Snyder,  and  S.  P.  Jacobs,  "Computational 
Considerations  for  Maximum-Likelihood  Radar  Imaging,"  Proc.  1990  CISS,  Princeton  Uni¬ 
versity. 


-  51  - 


Computational  Considerations  for  Maximum-Likelihood  Radar 

Imaging^ 

J.  A.  O'Sullivan 
P.  Moulin 
D.  L.  Snyder 
S.  P.  Jacobs 

Electronic  Systems  and  Signals  Research  Laboratory 
Department  of  Electrical  Engineering 
Washington  University 
Saint  Louis,  MO  63130 


ABSTRACT 

Recent  papers  have  outlined  a  new  approach  for  spectrum  estimation  and  radar 
imaging  based  on  expectation-maximization  algorithms  for  structured  covariance  estima¬ 
tion.  Performance  of  this  approach  has  been  promising  for  the  problems  studied.  Appli¬ 
cation  of  this  approach  to  real  data  sets  has  been  limited,  however,  due  to  the  need  to 
invert  a  matrix  whose  dimension  equals  the  size  of  the  data  set.  For  radar  applications 
where  an  image  is  to  be  formed,  data  sets  can  be  on  the  order  of  214  for  128x128  images. 

This  makes  the  use  of  the  new  approach  difficult  in  its  previously  described  form.  This 
paper  proposes  both  approximation  methods  for  inverting  typical  matrices  and  constraints 
on  radar  transmitted  signals  which  make  maximum  likelihood  image  estimation  viable. 

These  constraints  may  be  satisfied  for  real  signals  used  in  radar  imaging  systems.  Simu¬ 
lations  are  shown  to  demonstrate  the  performance  of  the  algorithms.  Finally,  motivated 
by  the  images  resulting  from  the  simulations,  regularization  methods  are  discussed. 

Introduction 

New  approaches  are  being  studied  for  maximum  likelihood  spectrum  estimation  and  radar  imaging 
which  arc  based  on  using  the  expectation-maximization  algorithm  [1]  for  structured  covariance  estimation 
[2-4).  In  this  paper,  we  focus  on  the  radar  imaging  problem,  although  many  of  the  results  hold  for  similar 
spectrum  estimation  problems.  The  limitations  of  our  previous  algorithm  [4]  arc  very  clear  for  the  radar 
imaging  problem  with  large  data  sets.  In  order  to  form  an  image  from  a  data  set  of  size  /V,  a  matrix  of 
dimension  NxN  must  be  inverted.  When  N  is  on  the  order  of  214,  this  inversion  is  not  practical  without 
exploiting  its  special  structure  or  making  some  approximations. 

First,  the  equations  which  describe  the  radar  data  arc  defined.  Next,  the  algorithm  derived  in  [4]  is 
presented.  After  discussing  the  role  of  the  matrix  inverse  in  the  algorithm,  possible  implementations  on 
massively  parallel  machines  arc  proposed.  Even  the  huge  number  of  processors  available  on  massively 
parallel  machines  cannot  make  the  inversion  problem  tractable  without  additional  assumptions.  For  practi¬ 
cal  radar  imaging  problems,  the  matrix  to  be  inverted  is  Toeplitz-Block  Toeplitz,  so  some  savings  in  com¬ 
putations  are  possible  by  exploiting  this  structure.  Other  improvements  are  possible  by  making  further 
assumptions  about  the  transmitted  signal  and  the  image  to  be  formed.  Finally  simulations  arc  presented 
and  the  need  for  regularization  methods  discussed. 

A  model  for  the  received  signal  in  a  radar  system  is  given  in  [5]  for  reflections  at  microwave  fre¬ 
quencies  and  in  (6)  for  reflections  at  optical  frequencies.  The  reflectivity  process  which  characterizes  the 
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target  is  a  random  process  which  is  stationary  at  each  delay.  Thus  the  samples  of  the  reflectivity  process 
have  Toeplitz  covariances,  and  the  problem  of  estimating  the  parameters  of  the  underlying  spectra  reduces 
to  a  Toeplitz  covariance  estimation  problem.  We  take  circulant  extensions  of  the  covariance  matrices  and 
estimate  spectrum  samples. 


1.  Statistical  Model 

The  target  is  described  in  terms  of  its  reflectivity.  It  is  assumed  that  the  radar  transmitted  signal  is  a 
plane  wave  at  the  target  so  points  on  the  target  at  a  cross-section  perpendicular  to  the  line  of  sight  sum  up 
to  contribute  to  the  same  return  signal.  This  sum  changes  as  a  function  of  time  and  is  denoted  by  b  ( t  ,x). 
The  variable  x  is  the  distance  of  points  on  the  target  from  the  transmitter  given  in  time  units  as  the  time  it 
takes  for  a  wave  to  propagate  to  the  target  and  back  to  the  transmitter  (two-way  delay).  From  this  model,  it 
is  not  apparent  that  separate  points  at  the  same  two-way  delay  x  may  be  differentiated  to  obtain  an  image. 
These  points  may  be  differentiated  if  their  velocities  relative  to  the  transmitter  are  different.  In  particular, 
if  the  target  is  a  rotating  rigid  body,  then  the  velocity  of  a  point  in  the  direction  of  the  transmitter  is  propor¬ 
tional  to  the  distance  of  that  point  from  the  line  of  sight.  Since  the  Doppler  shift  introduced  by  a  point  on 
the  target  is  proportional  to  this  directed  velocity,  a  delay-Doppler  image  of  a  rotating  rigid  body  is 
equivalent  to  a  range-crossrange  image.  The  problem  is  to  determine  the  power  reflected  from  points  as  a 
function  of  delay  and  Doppler  and  then  to  display  this  power  function  as  an  image  of  the  target.  Even  if 
this  rigid  body  assumption  for  the  target  is  not  valid,  a  delay-Doppler  image  can  be  a  useful  image  of  the 
target  region. 

The  reflectivity  b  (t  ,x)  is  the  superposition  of  all  reflectivities  at  delay  x,  times  the  Doppler  shift  terms 
they  introduce.  It  may  be  expressed  as 

b(l,x)  =  jc(f,x)ej2K/‘df,  (1.1) 

where  c  (f  ,x)  is  the  reflectivity  of  the  points  on  the  target  at  two-way  delay  x  which  introduce  a  Doppler 
shift  / .  The  target  is  assumed  to  have  finite  extent.  This  implies  that  both  b  ( t  ,x)  and  c  ( f  ,x)  are  zero  for  x 
outside  of  some  fixed  interval.  It  also  implies  that  c  (f  ,x)  is  zero  for  /  outside  of  some  finite  interval 
because  the  Doppler  variable  corresponds  to  crossrange  extent  of  the  target. 

The  literature  on  radar  reflections  [5,6]  describes  statistical  models  for  the  reflectivity  when  the  target 
is  rough  in  the  sense  that  multiple  scattering  sites  are  present  in  a  resolution  patch  on  the  target.  The  model 
states  that  the  reflectivity  of  a  patch  on  the  target  is  a  complex  valued  Gaussian  random  variable  with  zero 
mean  and  is  uncorrelated  from  patch  to  patch.  When  our  model  is  discretized  with  Icr^r  resolution  cells 
(the  subscripts  CR  and  R  denote  cross-range  and  range,  respectively),  this  corresponds  to  the  assump¬ 
tions  on  the  reflectivity  of  patches  of  the  target 

E[c(k,i)c{k\i'))  =  0 

E  [c  (*  ,i )  c *  (k  0]  =  a(k  J )  5U .  61( ;  • ,  (1.2) 

for  -(lCR-\)/2<k  <(ICr- l)/2,  0  <t  <!r-  1  .  Here,  a(k,i)  is  the  real,  nonnegative  covariance  of  the 
reflectivity  of  the  patch  at  Doppler  k  and  delay  /  .  The  scatlcrcrs  described  by  this  model  arc  called  widc- 
sense  stationary  uncorrelated  scattcrcrs  (WSSUS)  by  Van  Trees  [5]  because  the  assumptions  imply  that 

E  [b  (n  ,i  )b *  (n  ')]  =  Kb  (n  -n  \i  )5„  • ,  (1.3) 

where  b(n,i)  =  b  (n  At  -i  Ax/2  Ax)  and 

<Ja  - 1  >  2 

Kb(nj)=  2  c(k,i)cxp[j2KknfDAt\  ,  (1.4) 

*=~</o.-l)/2 

where  f  D  is  the  step  between  successive  Doppler  samples  in  the  image  plane.  Ax  is  the  delay  step,  and  At 
is  the  time  between  successive  samples  of  r(t). 

In  a  radar  system,  each  signal  is  typically  described  as  the  product  of  a  baseband  signal  times  a  com¬ 
plex  exponential  at  the  carrier  frequency.  All  of  the  interactions  of  interest  for  narrow  band  radar  systems 
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may  be  described  in  terms  of  these  complex  valued  baseband  signals  which  will  be  called  complex 
envelopes  of  the  signals  or  simply  the  signals. 

Let  sT(t)  and  sR(t)  be  the  complex  envelopes  of  the  transmitted  and  received  signals,  respectively. 
The  received  signal  is  the  superposition  of  the  reflectivity  from  all  two-way  delays  times  the  appropriate 
transmitted  signals, 

s„(t)=  ]sT«-x)b(t-x/2,x)dx.  (1.5) 

In  Equation  (1.5),  the  x/2  arises  because  it  takes  that  long  for  the  reflected  signal  to  return  to  the  receiver. 
The  available  data  are  the  sum  of  the  radar  return  signal  and  additive  white  Gaussian  noise 

r(t)~sR(t)  +  w(t) 

r  (r )  =  j  sT(t  -x)b  (t  -x/2,x)d  t  +  w  (r )  .  (1.6) 

This  equation  forms  the  basis  for  our  model. 

The  processing  is  assumed  to  be  performed  digitally  so  that  discretized  versions  of  b  and  c  are  used. 
Equation  (1.1)  is  substituted  into  (1.6)  and  approximated  by  samples  of  the  signals  c  and  w.  The  discrete 
equation  may  be  written  in  vector  form  as 

r  =  T^c  +  w,  (1.7) 

where  T  is  a  matrix  whose  entries  are  samples  of  the  transmitted  signal  times  appropriate  complex 
exponentials,  and  the  vectors  c  and  w  are  defined  appropriately. 


2.  Maximum-Likelihood  Solution 

The  loglikelihood  function  for  the  data  is 

L  (K*  ;r)  =  -ln(detK* )  -  rf{KR  )-'r,  (2.1) 

where  KR  is  the  covariance  matrix  for  the  received  data  r, 

k*  =  r'xr  +  jVqI,  (2.2) 

and  I  is  a  matrix  with  diagonal  entires  a(i,k).  As  shown  first  by  Burg,  et  al.  [3],  a  necessary  condition  for 
a  matrix  to  maximize  the  loglikelihood  function  is  the  trace  condition 

tr  [(KrWk*-1  -  K*')SKR  ]  =  0.  (2.3) 

The  matrix  8KR  is  a  variational  matrix  which  takes  values  in  all  possible  additive  variations  of  the  matrix 
Kfl ,  and  may  be  rewritten  as  T^ZT,  where  5Z  is  a  diagonal  variational  matrix. 

An  EM  algorithm  has  been  derived  [4]  to  estimate  the  matrix  E.  The  estimate  at  step  p+1  of  I  is 
given  by  the  diagonal  elements  of  the  conditional  expectation  of  ccf  or 

d 

£(p+i)  _  £  [ccf|Z(p),r],  (2.4) 

d 

where  =  means  that  the  off  diagonal  terms  on  the  left  are  zero  and  the  diagonal  terms  on  the  left  equal  the 
diagonal  terms  on  the  right.  The  computation  indicated  in  (2.4)  is  a  standard  problem  in  estimation  theory. 
This  equation  may  be  written  as  [4] 

d 

I^+1)  =  l'o)+  (2.5) 

i^rxr'i^’r  +  WoirV'  -  r'z^r  -  A/oixr'i^r  +  Woir'r'i^. 

Equation  (2.5)  defines  the  iteration  sequence  used  by  the  computations.  To  start  the  algorithm,  l{Cl)  is 
chosen  as  an  arbitrary  positive  definite  diagonal  matrix. 
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There  are  two  possible  images  which  may  be  displayed.  The  first  consists  of  estimates  of  samples  of 
the  scattering  function.  The  second  consists  of  magnitudes  or  squared  magnitudes  of  estimates  of  samples 
of  the  reflectivity  function.  The  diagonal  elements  of  I  are  the  values  which  are  displayed  as  the  scattering 
function  image  of  the  target.  Thus,  at  each  stage  an  image  is  calculated  and  may  be  displayed.  Some 
appropriate  stopping  criterion  is  used  to  terminate  the  algorithm.  At  each  stage  of  the  algorithm  the  condi¬ 
tional  mean  estimate  of  the  reflectance  is  also  generated.  At  step  p,  this  estimate  is 

£  Cellar]  =  j^)r(r,Z<p>r  +  A/0ir1r.  (2.6) 

The  magnitude  or  the  magnitude  squared  of  c  may  also  be  displayed  at  each  stage  of  the  algorithm  as  an 
image  of  the  target.  Thus  both  types  of  radar  image  commonly  viewed  are  generated  by  our  algorithm. 
We  feel  this  is  a  unique  feature  of  our  algorithm. 

Let  the  corresponding  sequence  of  covariance  matrices  for  the  data  r  be  denoted 
+  N0l.  Since  this  iteration  is  an  EM  algorithm,  it  has  all  of  the  properties  associated  with 
this  type  of  algorithm.  In  particular,  the  incomplete  data  loglikelihood  is  nondecreasing  in  the  sequence  of 
covariance  matrices  K jf\ 

There  are  obviously  issues  associated  with  the  appropriate  or  desired  sampling  rates  of  r(t)  and  of 
c(f , t).  Some  of  these  issues  are  addressed  in  [4],  The  quality  of  the  image  obtained  and  the  resolution 
achievable  are  intimately  related  to  the  sampling  issues. 

This  section  has  presented  a  review  of  the  equations  used  to  produce  target  images.  This  approach 
starts  with  a  model  which  accurately  accounts  for  the  random  nature  of  radar  reflections  and  adopts  the 
maximum  likelihood  method  of  statistics  to  estimate  delay-Doppler  high  resolution  images  of  radar  targets. 
Questions  associated  with  uniqueness  of  spectrum  estimates  and  convergence  of  the  EM  algorithm  arc  dis¬ 
cussed  in  [4,7],  Equation  (2.5)  is  a  computationally  demanding  update.  Some  of  the  issues  associated  with 
this  update  are  addressed  in  the  next  section. 


3.  Implementation  of  Algorithm 

The  implementation  of  the  iteradve  algorithm  described  in  (2.5)  involves  a  number  of  linear  algebra 
operations  on  arrays  which  may  be  very  large.  As  a  result,  the  development  of  efficient  and  numerically 
stable  routines,  on  both  serial  and  massively-parallel  machines,  to  perform  these  operations  is  of  great 
importance  in  evaluating  this  algorithm.  Because  massively-parallel  architectures  are  relatively  new  and 
higher  level  software  is  not  yet  available,  algorithmic  development  for  these  machines  is  more  involved. 

In  order  to  understand  the  issues  associated  with  our  effort  to  apply  the  algorithm  of  (2.5)  to  a  paral¬ 
lel  architecture,  one  must  understand,  to  a  degree,  the  limitations  imposed  by  our  resident  parallel  machine, 
the  DAP  510,  and  its  programming  language,  Fortran-Plus.  Fortran-Plus  is  an  adapted  version  of  Fortran 
which  hides  all  communication  between  processor  elements  from  the  user  by  virtue  of  its  data  structures. 
Scalars,  length-32  vectors,  and  32x32  matrices  are  all  distinct  data  types,  and  arrays  of  each  may  be 
formed.  Additionally,  Fortran-Plus  supports  no  complex  data  type.  Thus,  the  large,  complex  matrices 
found  in  our  signal  model  should  be  represented  as  3-dimensional  arrays  of  Fortran-Plus  matrices  in  order 
to  take  maximum  advantage  of  the  parallel  nature  of  the  DAP.  As  we  shall  see,  this  in  some  ways  compli¬ 
cates  and  in  other  ways  simplifies  the  task  of  programming. 

3.1  Algebraic  Operations  in  Algorithm 

If  one  examines  (2.5),  it  is  seen  that  the  indicated  iterations  can  be  performed  by  making  use  of  only 
the  following  basic  linear  algebra  operations: 


c  = 

A  +  B 

c 

=  AB 

c 

=  a'b 

A 

=  vwf 

w 

=  Afv 

-  55  - 


-5- 


B  =  A-1  (3.1) 

where  A,  B,  and  C  are  complex  matrices,  and  v  and  w  are  complex  vectors.  The  addition  operation  is 
greatly  simplified  by  the  data  structures  of  Fortran-Plus;  one  need  only  add  each  of  the  32x32  components 
of  A  to  the  corresponding  component  of  B.  Furthermore,  each  of  these  32x32  additions  can  be  accom¬ 
plished  by  a  single  statement  in  Fortran-Plus. 

Each  of  the  fear  multiplicadon  operations,  although  different,  can  be  performed  in  Fortran-Plus  by 
following  a  similar  strategy.  Consider  the  first  of  these,  that  of  performing  matrix  multiplication,  where  the 
matrices  in  question  are  complex  and  of  size  32/1x32/: ,  where  n  is  a  integer  greater  than  1.  Thus,  these 
matrices  would  be  represented  as  an  n  xn  x2  array  of  Fortran-Plus  matrices. 

Golub  and  Van  Loan  [8]  discuss  two  approaches  to  the  matrix  multiplication  problem  based  on  the 
following  partitioning: 

/zxnx2  (3.2) 

The  more  efficient  of  these  two  methods  is  that  developed  by  Strassen,  which  performs  the  multiplication 
of  two  M  xM  matrices  via  7  multiplies  and  18  adds  of  M  I2xM  /2  matrices: 

Pi  =  (An  +  A22XB11  +  B22) 

P2  =  (A21  +  Ayi)B]i 

P3  =  A„(Bi2  —  B22) 

P4  =  A22(B2i  -  Bu) 

P5  =  (An  +  A  12)822 

P6=  (A2I  -  An)(Bu  +  B12) 

P7  =  (Aj2  -  A22)(B21  +  B22) 

C„  =  Pi  +  P4  -  P5  +  P7 
Cl2  =  P3  +  P5 
^21  =  P2  +  P4 

C22=P1  +  P3-P2  +  P6  (3.3) 

This  method  is  particularly  useful  for  application  to  the  DAP,  if  used  recursively  to  reduce  the  problem  to 
that  of  multiplying  and  adding  a  series  of  32x32  real  matrices,  which  are  straightforward  and  efficient 
operations  in  Fortran-Plus.  A  slight  modification  of  Slrassen’s  algorithm  is  required  to  handle  the  separate 
real  and  imaginary  submatriccs,  but  the  recursive  partitioning  scheme  remains  intact. 

The  issue  of  matrix  inversion  within  the  algorithm  is  more  fundamental,  as  the  choice  of  algorithm  is 
not  obvious.  Press,  et.  al.  [16]  discuss  the  merits  of  Strassen  inversion,  which  is  based  on  an  analogous  par¬ 
titioning. 

32x32  (3.4) 

R,  =  Aj]1 

R2 =  A2iXR, 

R3 =  R  i  x  A 1 2 

R4  =  A2ixR3 

R5  =  R4  -  A 22 

r6  =  r5-' 

C12  =  83*85 
c21  =  86xR2 
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R?  -  R3XC21 

C22  =  -R6  (3.5) 

The  recursive  partitioning  associated  with  this  algorithm  makes  it  particularly  attractive  for  implementation 
on  the  DAP,  as  all  numerical  operations  can  be  performed  on  Fortran-Plus  matrices. 

3.2  Further  matrix  inversion  issues 

In  addition  to  the  estimated  covariance  matrix  Kk  ,  often  the  matrix  T  is  inverted  for  the  purpose  of 
comparing  the  results  of  the  EM  algorithm  with  those  obtained  from  (1.7)  in  an  attempt  to  reconstruct  c 
using  linear  processing.  This  is  in  the  spirit  of  conventional  processing  which  is  based  on  linear  operations. 
In  addition,  this  corresponds  to  the  maximum  likelihood  solution  shown  in  (4.6)  in  the  special  case  dis¬ 
cussed  there.  Note  that  the  parallel  algorithm  for  inverting  matrices  in  (3.5)  seeks  to  invert  the  upper  left 
block  of  the  original  matrix.  It  happens  that,  under  certain  conditions,  the  upper  left  block  of  the  matrix  T 
may  be  singular,  or  at  least  poorly  conditioned.  As  the  partitioning  is  repeated,  the  algorithm  becomes 
numerically  unstable  for  a  large  class  of  matrices  T.  Other  algorithms  which  first  partition  the  matrix  to  be 
inverted  into  its  real  and  imaginary  parts  suffer  the  same  fate. 

If  one  examines  the  literature,  one  finds  a  wide  variety  of  inversion  algorithms  based  on  decomposi¬ 
tions.  Each  seeks  to  make  the  problem  more  tractable  by  factoring  the  matrix  in  question  in  such  a  way 
that  the  factors  are  easily  inverted.  Many  of  these  decompositions  require  that  the  matrix  to  be  inverted 
have  a  special  structure  which  either  Ks  or  T  do  not  satisfy.  For  example,  the  LU  factorization  requires 
that  all  principal  submatrices  of  the  matrix  in  question  be  non-singular,  in  order  for  the  algorithm  to  be 
stable  [8].  As  we  have  already  seen,  this  may  not  be  true  for  T. 

We  have  chosen  to  implement  the  QR  decomposition  in  Fortran-Plus  for  this  purpose.  This  routine 
factors  a  square  matrix  A  in  in  the  following  way: 

A  =  QR  (3.6) 

such  that  Q  is  unitary  (QfQ  =  I)  and  R  is  upper  triangular.  This  gives  rise  to  the  following  inversion  algo¬ 
rithm: 

A-1  =  R_1Qf  (3.7) 

Thus  the  inversion  problem  has  been  transformed  to  the  separate  problems  of  calculating  the  factors, 
inverting  a  triangular  matrix,  and  performing  one  matrix  multiplication.  Note  that  although  this  process 
docs  not  reduce  the  complexity  of  the  inversion  problem,  it  defines  the  solution  in  terms  of  operations 
which  are  well  known. 

Unfortunately,  the  parallel  implementation  of  this  inversion  strategy  is  not  as  computationally 
efficient  as  one  might  hope,  as  the  operations  associated  with  the  QR  decomposition  can  really  only  be 
applied  to  columns  of  data  at  a  time.  In  order  to  take  maximum  advantage  of  the  DAP’s  parallel  architec¬ 
ture,  such  operations  should  be  applied  to  32x32  matrices  at  a  time.  However,  the  inverse  of  T  is  computed 
only  once  for  the  purposes  of  forming  the  output  of  conventional  processing,  while  the  inverse  of  is 
computed  at  every  iteration  of  the  EM  algorithm.  Therefore,  a  significant  increase  in  the  speed  of  the  algo¬ 
rithm  as  a  whole  may  be  realized  if  a  more  efficient  parallel  routine  for  computing  is  utilized.  Certain 
special  cases  for  the  transmitted  signal  give  rise  to  such  situations,  as  discussed  in  the  next  section. 


4.  Improvements  in  Complexity 

From  the  considerations  developed  in  Section  3,  it  appears  that  the  task  of  inverting  the  covariance 
matrix  K*  is  a  formidable  one.  In  this  section  we  examine  under  what  conditions  the  complexity  of  this 
operation  is  lower  than  first  expected. 

4.1  Toeplitz-Block  Toeplitz  structure 

Consider  the  often-used  steppcd-frequency  waveform  [9].  The  transmitted  signal  is  made  of  Nb 
bursts  of  Np  pulses  each.  The  signal  transmitted  in  pulse  p„  of  burst  b„  is  a  complex  sinusoid  at  frequency 
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pJT  (where  fT  is  the  frequency  step  between  successive  pulses): 

sT(0  =  exp[-jlxpJTt]  ,  (4.1) 

Usually  the  return  signal  is  sampled  at  the  pulse  repetition  frequency  (1  sample  per  pulse  ;  N  =  NbNp).  If 
f  rAt  is  an  integer,  we  obtain 

sT(n  At  -i  Ax)  =  exp {j  2 npnf  Ti  Ax] ,  (4.2) 

with  n  =  Npbn  +  pn,  0<bn  <Nb  ,0<pH  <Np.  In  this  instance  K*  has  a  special  structure,  known  as 
Toeplitz-Block-Toeplitz  (TBT)  [10].  From  (2.2): 

/d»  — 1/»— 1 

K R(n,m)=  £  X •  * ) expty 2nkfD (n -m )At ] sT(n At-i Ax) sj(m At -i Ax) 

k=0  i=Q 

Zo.-1/.-l 

=  £  £  o(k ,  i )  exp  [j  2nkfD  (n  -m  )A t  ]  exp  [j  2  n(pH  -pm  )fTi  Ax] 

*=0  i=0 

=  X  Tia(.k’i)<ixp[j2nkfD(Np(bn-bm)+(plt-pm))At]x 

i=0  i=0 

exp[j  2n(p„ -pm  )fTi  Ax] .  (4.3) 

Thus  K R(n,m)  is  a  function  of  bn-bm  and  p„-pm  only.  possesses  a  doubly  Toeplitz  structure:  Ks  is 
made  of  Nb  xNb  Toeplitz  blocks  of  dimension  Np  xNp ,  themselves  arranged  in  a  Toeplitz  structure. 

Special  algorithms  have  been  developed  for  inverting  TBT  matrices  [10,11],  The  complexity  of 
Wax’s  algorithm  is  min  (NpNb,  NpNb).  For  Nb  =Np  this  is  equal  to  N5'2,  which  is  significantly  smaller 
than  N3  for  large  N .  The  complexity  of  the  TBT  matrix  inversion  problem  on  a  parallel  machine  has  not 
yet  been  investigated. 

4.2  Special  choice  of  the  parameters 

In  the  special  case  where  IrIcr  =N  and  ITf  is  the  identity  matrix,  the  covariance  matrix  (2.2)  has 
the  form 

KR=rt(l+N0l)r.  (4.4) 

The  trace  condition  (2.3)  becomes 

7V  [5m+W0ir1(yyf-^V0I)(I+N0ir1]  =  0  ,  (4.5) 

where 

y  =  (r')fr  .  (4.6) 

The  covariance  to  be  inverted  in  this  equation  is  diagonal.  Numerically  the  inversion  problem  is  now 
trivial.  Actually  this  equation  can  even  be  solved  in  closed  form.  The  maximum-likelihood  estimates  are 
given  by: 

a(k,i)  =  max  |y  (k,i)\2-N0  ,  oj  .  (4.7) 

It  should  be  noted  that  for  stepped- frequency  waveforms,  the  design  ICR  =  Nb  and  1R  =  Np  leads  to  rrf 
being  the  identity  matrix.  In  this  instance,  |y  |2  is  the  output  of  the  conventional  processing  for  the  radar 
return,  based  on  Fourier  transforms  [9],  The  ML  estimator  is  obtained  by  subtracting  ,V0  from  this  estimate, 
followed  by  setting  negative  values  to  zero. 

5.  Regularization  Methods 

The  ML  estimator  presented  in  Section  2  has  been  shown  to  offer  lower  bias  and  better  resolution 
than  conventional  estimators.  However  the  estimates  obtained  with  the  estimation  setup  of  Section  2  arc 
not  mean-square  consistent.  This  unfortunate  result  is  due  to  the  the  number  of  parameters  (  IcrIr  )  being  at 
least  equal  to  the  number  of  data  ( N )  in  the  method  used.  We  arc  actually  dealing  with  a  small-sample 
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MLE  problem.  The  situation  becomes  worse  as  IrIcr  is  increased  in  an  attempt  to  improve  resolution. 
The  estimates  are  extremely  sensitive  to  small  changes  in  the  data  and  exhibit  a  very  rough  shape  charac¬ 
terized  by  very  sharp  peaks  and  low  valleys. 

Such  a  behavior  is  of  course  highly  undesirable,  and  we  need  a  regularization  of  the  estimates.  This 
concept  was  introduced  by  Tikhonov  [12]. 

One  possible  approach  for  regularization  of  the  estimates  is  the  method  of  sieves  introduced  by 
Grenander  [13,  p.  357].  A  sieve  in  a  parameter  space  A  is  a  family  of  subsets  S(p)  of  A  indexed  by  a  posi¬ 
tive  parameter  p.  called  the  mesh  size.  A  restricted  ML  estimate  exists  over  each  set  S(p).  As  the  mesh  size 
tends  to  zero  the  sets  S(p)  will  be  large  enough  to  allow  the  ML  solution  to  converge  to  any  solution  in  A. 
For  the  problem  at  hand  the  choice  of  a  "good"  value  of  p  will  depend  on  the  data  record  size  and  the  noise 
level,  among  other  factors. 

A  major  problem  is  to  find  "good"  sieves,  which  will  make  the  ML  estimates  converge  and  possess 
some  desirable  practical  features  such  as  analytical  and/or  computational  tractability.  Recently  we  pro¬ 
posed  a  method  of  sieves  for  a  spectrum  estimation  problem  [14],  The  function  to  be  estimated  is 
represented  by  a  series  expansion,  and  the  restricted  set  S(Q)  is  the  set  of  series  truncated  to  Q  terms. 
Essentially  the  variational  matrix  8 KR  is  now  constrained  in  a  Q  -dimensional  subset  which  will  be  much 
smaller  than  the  set  of  all  possible  variations.  This  method  can  be  extended  to  the  radar  imaging  problem 
as  follows  [15]. 

5.1  Definition  of  a  Sieve 


The  scattering  function  a(k,i)  can  be  viewed  as  a  vector  in  the  Euclidean  space  R/c,/’.  Let 
y/m  {k ,  i ) ,  0  <  /  <ICr  ,  0  <  m  <  IR  ,  be  a  set  of  ICrIr  positive  basis  vectors  in  .  It  is  not  required 
that  this  set  be  orthogonal.  Consider  the  set  of  all  vectors  in  the  span  of  ]  with  nonnegative 
coefficients: 


/«- 1  /.- 1  l 

A  :=  <  <5{k ,  i )  =  2  E  a(l,m)\Vim^’i)  ,a(l,m)>  Ok 
I  1=0  m= 0  J 


(5.1) 


We  define  our  parameter  set  to  be  A.  The  vectors  in  this  set  are  represented  by  means  of  the  coefficients 
a(l,m)  0<l  <1Cr  ,  0<m  <IR.  We  define  a  sieve  S (Qcr<Qr)  on  a (k,i)  by  truncating  the  scries 
representation  (5.1)  to  QCrQr  terms.  The  integers  QCR  and  QR  are  constrained  to  be  no  larger  than  lCR 
and  Ir,  respectively. 

f  Qcn-l  Qm-1  1 

S (QCR’QrY-=  \  a(k,i)=  2  'L  a(l<m)ylm(k,i)  ,a(l,m)>0>  .  (5.2) 

/  =^J  m= 0  J 

The  vectors  in  this  subset  are  represented  by  means  of  QcrQr  coefficients 
a  (l  ,m)  0<l<  Qcr  ,  0  <m  <QR.  We  view  these  coefficients  as  being  QCR  QR  unknown  parameters  for 
which  ML  estimates  are  sought.  In  [15]  we  show  that  consistency  of  the  estimates  can  be  obtained  if  QCr 
and  Qr  grow  at  an  appropriate  rate  with  N .  Typically  the  basis  functions  that  are  designed  arc  localized  in 
range  and  cross-range  (i.e.  k  and  i).  Next  we  show  how  the  ML  estimates  of  the  coefficients  arc  com¬ 
puted. 


5.2  Algorithm  No.l 

For  each  l,m ,  define  the  support  set  D(m  of  \\f^,(k,i).  Also  construct  a  diagonal  matrix  4/,m  made  of 
the  samples  i|/;m (k,i)  and  a  basis  covariance  matrix 

K,m:=rf%mr.  (5.3) 

It  follows  that  Kr  ,  subject  to  the  sieve  constraint,  assumes  the  form 

Qc-\Qa~\ 

K*  =  I  I  a(l,m)Kln  +N0 1.  (5.4) 

1=0  m  =0 

The  estimation  problem  is  now  to  maximize  the  loglikelihood  (2.1)  subject  to  the  constraint  (5.4).  From 
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(2.3)  and  (5.4),  a  necessary  condition  for  Ks  to  maximize  (2.1)  is  the  trace  condition 

tr  [  (K« ‘rr'K*1  -  K*1)  Klm  ]  =  0  ,  0  <  l  <  QCr  ,  0  <  m  <  Qr  .  (5.5) 

The  maximization  problem  is  solved  using  an  EM  algorithm  similar  to  the  one  derived  in  [14],  The  com¬ 
plete  data  are  defined  as  ({ctm},  w  ),  where  clm  ,  0</  <  QCR  ,  0<m  <QR,  is  an  ICrIr -vector  with 
covariance  a(/,m)vP/m  known  up  to  the  scale  factor  a(l,m).  The  estimate  at  step  p+1  of  the  coefficient 
a  (/ ,  m )  is  given  by 

w„+n  1  ^  £[|Ch,(*.i)|2k.<f(/.i«)w] 


d(/,/n)(p+1)  = 


I  k,ie  Dta 


Evaluadng  the  expectation  yields  the  update  equation 


T~  [  <f(/.«H  2  tr  [(K^^rr^^-K*1)^] . 

I  Vlm  I 


The  algorithm  is  started  with  positive  coefficients  a(l,m).  Each  matrix  K,m  is  constant  and  can  be  com¬ 
puted  off-line.  At  each  step  the  covariance  matrix  Kjf)  must  be  inverted,  which  accounts  for  JV3  opera¬ 
tions.  Following  this,  updating  each  a(l,m)  just  requires  N2  multiplications/additions.  The  complexity  of 
each  step  of  the  EM  algorithm  is  then 

N3  +  QcrQrN2  ,  (5.8) 

as  compared  to  A3  +  IcrIrN1  in  the  unconstrained  case  of  §2.  This  saving  is  important,  as  it  is  now  possi¬ 
ble  to  let  ICR IR  tend  to  infinity.  However  the  problem  of  inverting  the  covariance  matrix  is  just  as  difficult 
as  before.  In  the  next  section  we  show  how  the  method  is  handled  in  a  special  case  of  interest. 

5.3  Algorithm  No.2 

The  regularization  method  described  in  §5.2  can  be  applied  to  the  special  case  presented  in  §4.2.  In 
general  the  trace  condidon  cannot  be  solved  in  closed  form,  yet  the  computauonal  requirements  are  drasu- 
cally  reduced.  The  update  equadon  (5.7)  becomes 

d(l,m)(p*l)  =  d(l,m)<J,)+  — [  a(/,m)(p)  1  2x 

I  D/m  I  ^  J 

tr  [((I(P)  +  W0I)-1yy^)  +  /V0ir1-(Z^)  +  A0I)-1)^  j 

If  2 

=  d(/,m)(p)+  — ! —  x  (5.9) 

I  D/m  I  L 

\y(k,i)\2-a(k,i)^-N0 

,  ..  .  .  2  Vlm(k  ,l)  . 

k.ieD*  (G(k,l)V>  +  N0)2 

where  the  vector  y  is  defined  in  (4.6).  The  matrix  KR  has  been  reduced  to  a  diagonal  form  and  the  numeri¬ 
cal  complexity  of  each  iteration  is  equal  to 

Qa, —  1  Qjt  — 1 

I  IlD/m! 

/= 0  m= 0 

muldplications/additions.  This  number  is  the  total  area  of  the  support  sets  of  the  basis  funcuons.  For  a  dis¬ 
cussion  on  the  design  of  the  basis  funcuons  and  their  support  sets  we  refer  to  [15].  Typically  it  is  required 
that  the  support  sets  cover  the  entire  range  and  cross-range  domain,  so  that  the  complexity  is  equal  to 

Gc.-ifi.-i 

I  I|D,ml  =  McrIr  =  M  .  (510) 

/  =0  m- 0 

where  the  constant  X.  >  1  indicates  the  amount  of  overlapping  of  the  support  sets.  This  figure  is  a  major 
improvement  from  the  corresponding  figure  in  (5.8). 
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The  special  case  that  was  treated  in  this  section  is  very  important,  for  it  can  be  shown  that  the  esti¬ 
mates  obtained  with  the  algorithms  of  Sections  5.2  and  5.3  using  different  values  for  IcrIr  are  asymptoti¬ 
cally  equivalent  as  N  — > «  [15],  This  result  justifies  the  choice  of  setting  IcrIr  =  N ,  followed  by  using  the 
latter  algorithm,  as  a  valid  approach  to  estimating  a(k,i). 

6.  Simulations 

We  have  performed  simulations  which  will  allow  us  to  compare  the  results  of  the  algorithm 
described  in  section  2,  the  algorithm  of  section  5,  as  well  as  the  results  of  conventional  processing.  In  the 
first  of  these,  we  designed  a  series  of  routines  for  a  massively  parallel  machine  which  would  generate  simu¬ 
lated  radar  data  according  to  the  indicated  diffuse  target  model,  and  then  implement  the  algorithms  of  sec¬ 
tion  2  to  estimate  the  target  image.  Our  results  here  were  quite  promising.  In  every  case  examined,  our 
algorithm  outperformed  the  conventional  processing,  and  was  especially  superior  in  those  cases  where  the 
noise  level  was  significant. 

There  are  two  sets  of  plots  shown  below.  They  are  computed  for  a  scattering  function  which  is  a 
square  of  four  by  four  pixels  centered  in  a  square  of  size  eight  by  eight  pixels.  The  intensities  of  the 
scattering  function  in  each  of  the  regions  and  of  the  noise  are  the  parameters  in  the  plots  shown.  The 
transmitted  signal  is  a  stepped  frequency  waveform,  with  one  sample  taken  per  pulse.  A  total  of  64  pulses 
are  transmitted  in  the  simulations  with  eight  bursts  of  eight  pulses  each.  The  frequency  step  is  chosen  such 
that  fTA x  equals  0.125. 

The  first  plots  show  the  convergence  of  the  value  of  the  loglikelihood  function  for  single  realizations 
of  the  simulation  at  two  background  noise  levels.  The  signal  to  noise  ratio  indicated  on  the  plots  is  com¬ 
puted  by  first  taking  the  ratio  of  the  expected  total  energy  in  the  signal  part  of  the  received  signal  and  divid¬ 
ing  by  the  expected  total  energy  in  the  noise,  and  then  computing  10  times  the  log  of  that  number. 

The  second  set  of  plots  shows  the  summed  absolute  error  between  the  maximum  likelihood  image 
assuming  all  of  the  data  (c)  is  observed  and  the  maximum  likelihood  image  assuming  r  is  observed.  These 
plots  measure  the  degradation  in  performance  due  to  the  noise  and  having  to  invert  I~.  These  plots  show 
that  the  actual  values  displayed  in  the  image  continue  to  change  even  after  the  likelihood  function  is  close 
to  its  maximum.  Visually,  this  effect  is  manifested  in  the  image  becoming  rougher  as  the  iterations 
proceed.  The  regularization  methods  discussed  in  section  5  and  in  [15]  take  into  account  more  prior  infor¬ 
mation  known  about  the  scattering  function  and  thus  produce  less  rough  images. 
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HIGH  RESOLUTION  RADAR  IMAGING  USING  SPECTRUM  ESTIMATION  METHODS 


JOSEPH  A.  O'SULLIVAN 
DONALD  L.  SNYDER  t 


Abstract  This  paper  summarizes  a  new  approach  to  high  resolution  radar  imaging  based  on  modem 
spectrum  estimation  techniques.  First  a  statistical  model  of  the  radar  reflections  which  properly  accounts 
for  the  randomness  of  reflections  by  targets  which  are  rough  on  the  order  of  a  wavelength  of  the  carrier 
frequency  is  introduced.  The  model  for  the  radar  return  signal  is  valid  for  all  transmitted  narrowband  radar 
signals.  Equations  which  generate  maximum  likelihood  estimates  for  the  reflectivity  power  as  a  function  of 
delay  and  Doppler  coordinates  are  derived. 


Introduction.  This  paper  presents  recent  research  results  in  high  resolution  delay-Doppler  radar 
imaging  based  on  statistical  models  obtained  from  the  underlying  physics  of  spotlight  mode  radar 
reflections.  We  present  solution  equations  which  may  be  used  to  process  radar  return  signals  to  form 
images.  These  equations  bear  a  strong  resemblance  to  the  equations  derived  in  [1]  for  the  generic  Toeplitz 
covariance  estimation  problem.  This  is  a  result  of  the  fact  that  the  reflectivity  process  which  characterizes 
the  target  is  a  random  process  which  is  stationary  at  each  delay.  Thus  the  samples  of  the  reflectivity 
process  have  Toeplitz  covariances  and  the  problem  of  estimating  the  parameters  of  the  underlying  spectra 
reduces  to  a  Toeplitz  covariance  estimation  problem.  We  take  circulant  extensions  of  the  covariance 
matrices  and  estimate  spectrum  samples.  There  are  several  aspects  of  this  problem  which  differ 
significantly  from  the  generic  problem,  but  the  underlying  problem  is  to  estimate  Toeplitz  covariance 
matrices. 

The  model  of  the  reflectivity  process  is  described  first.  Since  the  processing  is  performed  digitally, 
the  discrete  form  of  this  model  is  examined  in  detail.  Next,  the  manner  in  which  the  transmitted  signal 
interacts  with  the  reflectivity  to  form  the  radar  return  signal  is  presented.  The  imaging  problem  is  reduced 
to  the  problem  of  estimating  the  spectrum  underlying  the  reflectivity  process  given  samples  of  the  radar 
return  signal.  A  necessary  condition  for  the  maximum  likelihood  solution  is  obtained  and  an  EM  algorithm 
approach  to  solving  for  the  maximum  is  taken.  The  results  extend  those  in  [1]  in  three  important  ways. 
First,  the  samples  of  the  process  with  a  Toeplitz  covariance  are  not  available.  Instead,  the  stationary 
process  is  multiplied  by  a  signal  matrix.  Second,  the  model  includes  an  additive  white  Gaussian  noise. 
Third,  the  process  of  interest  is  a  function  of  two  variables.  For  each  value  of  the  delay  variable,  the 
process  is  Toeplitz.  Thus  the  model  from  [1]  is  extended  to  spectra  which  change  as  a  function  of  an 
independent  variable. 


Reflectivity  process.  The  target  is  described  in  terms  of  its  reflectivity.  It  is  assumed  that  the  radar 
transmitted  signal  is  a  plane  wave  at  the  target  so  points  on  the  target  at  a  cross-section  perpendicular  to  the 
line  of  sight  sum  up  to  contribute  to  the  same  return  signal.  This  sum  changes  as  a  function  of  time  and  is 
denoted  by  y  (r  ,t).  The  variable  x  is  the  distance  of  points  on  the  target  from  the  transmitter  given  in  time 
units  as  the  time  it  takes  for  a  wave  to  propagate  to  the  target  and  back  to  the  transmitter  (two-way  delay). 
From  this  model  it  is  not  apparent  that  separate  points  at  the  same  two-way  delay  x  may  be  differentiated  to 
obtain  an  image.  These  points  may  be  differentiated  if  their  velocities  relative  to  the  transmitter  are 
different  In  particular,  if  the  target  is  a  rigid  body  and  is  rotating  about  a  point  along  the  line  of  sight  then 
the  velocity  of  a  point  in  the  direction  of  the  transmitter  is  proportional  to  the  distance  of  that  point  from 
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the  line  of  sight  Since  the  Doppler  shift  introduced  by  a  point  on  the  target  is  proportional  to  this  directed 
velocity,  a  delay-Doppler  image  of  a  rotating  rigid  body  is  equivalent  to  a  range-crossrange  image.  The 
problem  is  to  determine  the  power  reflected  from  points  as  a  function  of  delay  and  Doppler  and  then  to 
display  this  power  function  as  an  image  of  the  target  Even  if  this  rigid  body  assumption  for  the  target  is 
not  valid,  a  delay-Doppler  image  can  be  a  useful  image  of  the  target  region. 

Since  the  reflectivity  y(/,t)  is  the  superposition  of  all  reflectivities  at  delay  x  times  the  Doppler  shift 
terms  they  introduce,  it  may  be  expressed  as 

=  (1) 

where  c{f  ,x)  is  the  reflectivity  of  the  points  on  the  target  at  two-way  delay  t  which  introduce  a  Doppler 
shift  / .  The  target  is  assumed  to  have  finite  extent.  This  implies  that  both  y  (r  ,x)  and  c  (f  ,x)  are  zero  for  x 
outside  of  some  fixed  interval.  It  also  implies  that  c(f , t)  is  zero  for  /  outside  of  some  finite  interval 
because  the  Doppler  variable  corresponds  to  crossrange  extent  of  the  target  The  processing  is  assumed  to 
be  performed  digitally  so  that  discretized  versions  of  y  and  c  are  used.  Suppose  that  the  resolution  cells  of 
/  and  x  are  A/  and  At  respectively.  Let  there  be  /*  bins  in  the  delay  or  range  direction  (in  delay 
coordinates,  the  target  is  of  length  Ir  Ax/2)  and  let  there  be  Icr  bins  in  the  Doppler  or  cross  range  direction 
(so  the  target  is  of  width  lot  A/).  If  samples  of  the  radar  return  signal  are  taken  every  At  seconds. 
Equation  (1)  may  be  approximated  by 

(/«-  1V2 

y(kAtjiAx)  =  T  c(m,n)e/2*'"±V*  ,  (2) 

where 


c  (m  ji )  =  c  (m  A/  a  Ax)  A/.  (3) 

The  literature  on  radar  reflections  (2,3)  describes  statistical  models  for  the  reflectivity  when  the 
target  is  rough  on  the  order  of  a  wavelength  of  the  carrier  signal.  The  model  states  that  the  reflectivity  of  a 
patch  on  the  target  is  a  complex  valued  Gaussian  random  variable  with  zero  mean  and  is  uncorrelated  from 
patch  to  patch.  For  our  model  this  corresponds  to  the  assumptions  on  the  reflectivity  of  patches  of  the 
target 

E[c(i  Jc)c(mji)]  =  0 

E[c(ijc)c’(m/i)]  =c(iJx)5;jn5kA.  (4) 

Here,  a(i  Jc)  is  the  real,  nonnegative  covariance  of  the  reflectivity  of  the  patch  at  delay  k  and  Doppler  i . 
The  scatterers  described  by  this  model  are  called  wide-sense  stationary  uncorrelated  scatterers  (WSSUS) 
by  Van  Trees  [2]  because  the  assumptions  imply  that 

£1 j(M)y*(/«  J))  =  ^c(«-«^)8*j,  (5) 

where  y(iJc)  =  y(iAl  -k  Ax/2,Ax)  and 

KG(nJc)  =  Y  ■  (6) 

Let  yc  (i)  be  the  vector  of  samples  from  delay  k 

yc(*)'  =  [y«U)  y(U)  y{2Jc)  •••  y(G-U)].  (7) 

The  covariance  matrix  for  yo(k)  is  a  Toeplitz  matrix  Kc(k)  whose  ijn  element  is  The 

restriction  imposed  by  the  constraint  of  a  target  of  finite  extent  is  represented  by  Equation  (6).  That  is,  the 
only  a(mjc)  which  can  be  nonzero  are  for  -{Icg-l)/2Sm  S  (/c*-l)/2  and  0<k  £/#- 1.  This  constraint 
on  m  is  only  meaningful  if  A f  Allot  is  less  than  one  and  it  limits  the  Kc(jfc)  possible.  If  A/  A/  =  l/N ,  then 
the  a(mjc)  may  be  thought  of  as  samples  of  the  spectrum  of  the  periodic  extension  of  the  covariance 
matrix  at  delay  k .  Only  Icr  ( Icr  <N)  of  the  spectrum  samples  at  each  delay  are  nonzero. 

We  now  introduce  larger  vectors  and  matrices  so  that  the  results  which  follow  may  be  written 
compactly.  The  yc(k)  are  loaded  into  one  GIrxI  vector  y©.  The  covariance  of  yc  is  a  block  diagonal 
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matrix  Kg  whose  kA  diagonal  block  is  K a(k).  Define  I(k)  to  be  the  1Cr *Icr  diagonal  matrix 

Hk)  =  diag[a(0Jc)  a(U:)  o(2jk)  •••  a(-2Jc)  a(-ljc)  ]  (8) 

and  define  Z  to  be  the  IrIcr'xIrIcr  block  diagonal  matrix  whose  kA  block  is  Z(k).  Then  we  have 

Kg(*)  =  W£J&«*)Jc*Wg  (9) 


where  W  is  an  NxN  normalized  DFT  matrix  (N  >  G ,N  >  Icr),  Wc  consists  of  the  first  G  columns  of  W, 
and 


Jc* 


It 


0 


(10) 


where  Jc*  is  /c*xlV.  It  is  an  identity  matrix  of  dimension  (/c*+l)/2,  and  I2  is  an  identity  matrix  of 
dimension  (/c*-l)/2.  Jc*  implements  the  assumption  that  some  spectrum  samples  are  zero.  Having  now 
defined  several  matrices,  with  a  little  more  notation  we  can  give  a  more  precise  relationship  between  the 
matrix  Zand  the  matrix  Kg.  Let  Me*  beth  e/*/c*x/*iV  block  diagonal  matrix  each  of  whose  /*  blocks  is 
Jc*.  Let  Mr  be  the  IrNxIrG  block  diagonal  matrix  each  of  whose  /*  blocks  is  J*  (J*  is  NxG  and 
equals  [  I  0  ]')■  Finally,  let  W  be  the  NIrxNIr  block  diagonal  matrix  each  of  whose  /*  blocks  equals  W. 
Then 


Kg  =M*W'M£*ZMc*WM*. 


(11) 


Except  for  constraining  some  spectrum  samples  to  be  zero,  if  samples  of  yc(k)  from  each  delay  k 
were  directly  available,  this  problem  would  be  very  similar  to  that  from  [1]  and  the  approach  to  the  solution 
almost  identical.  Instead,  for  our  problem,  we  have  the  data  r  which  is  described  below. 


Radar  data.  In  a  radar  system,  each  signal  is  typically  described  as  the  product  of  a  baseband  signal 
times  a  complex  exponential  at  the  carrier  frequency.  All  of  the  interactions  of  interest  for  narrow  band 
radar  systems  may  be  described  in  terms  of  these  complex  valued  baseband  signals  which  will  be  called 
complex  envelopes  of  the  signals  or  simply  the  signals. 

Let  stU)  and  s*(r)  be  the  complex  envelopes  of  the  transmitted  and  received  signals,  respectively. 
The  distance  of  a  point  on  the  target  from  the  transmitter/receiver  is  measured  in  time  units  as  the  two-way 
delay.  The  received  signal  is  the  superposition  of  the  reflectivity  from  all  two-way  delays  times  the 
appropriate  transmitted  signals, 

Sr( t)-  JsT(t-x)y(t^i/2,x)dX.  (12) 

In  Equation  (12),  the  x/2  arises  because  it  takes  that  long  for  the  transmitted  signal  to  get  to  a  point  at  two- 
way  delay  x.  The  available  data  are  the  sum  of  the  radar  return  signal  and  additive  white  Gaussian  noise 

r(0  =  SR(t)  +  w(t) 

r(t)  =  J'ST(t^x)y(t-x/2,x)dx  +  w(t)  .  (13) 

This  equation  forms  the  basis  for  our  model.  We  now  assume  that  samples  of  the  data  are  available  and  we 
wish  to  estimate  the  covariance  of  the  reflectivities  from  patches  on  the  target.  The  sampled  version  of 
Equation  (13)  is 

r(k)  =  ^ST(kAt-nAx)y(k&t-nAx/2ji£x)  +  w(k)  .  (14) 

Given  G  samples  r  (k).  Oik <G - 1 ,  we  may  write  Equation  (14)  in  vector  form  as 

r  =  S^y  +  w  ,  (15) 

where  r  is  the  vector  of  samples  r(k);  w  is  a  white  noise  vector  with  covariance  Nol;  y  is  the  G/*x  1  vector 
of  samples  of  the  reflectivity  of  the  target  described  above;  and  Sf  is  a  G  *GIr  matrix  composed  of  /* 


-  68  - 


(16) 


G  xG  submatrices  each  of  which  is  diagonal: 

S*-[So  Si  Si  •••  S/,-1  ], 

where 

S„  =diag[ sri-nAx)  sriAt-nAx)  sr(2At-nAx)  •••  st((G -l)A:-n  Ax)  ]. 

Note  that  Equation  (15)  could  also  be  written  in  terms  of  samples  of  c  (/  ,t)  as 

r  =  Fc  +  w,  (17) 

where  c  is  an  /*/arX  1  vector  of  samples  of  the  reflectivity  of  the  target  arranged  as  Ir  subvectors  each  of 
length  lex  of  the  samples  from  each  delay  and  Tf  is  a  Gx/rIcr  matrix  each  element  of  which  is  a  product 
of  a  sample  of  st  times  a  complex  exponential.  More  explicitly,  T  is  given  by 

r  =  McxWM*Sr.  (18) 

There  are  obviously  issues  associated  with  the  appropriate  or  desired  sampling  rates  of  r(t)  and  of 
c(f  ,i).  Some  of  these  issues  are  addressed  in  [4].  The  quality  of  the  image  obtained  and  the  resolution 
achievable  are  intimately  related  to  the  sampling  issues. 

Since  r  is  the  result  of  linear  operations  on  Gaussian  random  variables,  r  is  a  zero  mean  Gaussian 
random  variable  with  covariance 

Kk  =SfKGST  +iVoI  =  rZT  +  NoI,  (19) 


Maximum  likelihood  solution.  The  loglikelihood  function  for  the  data  is 

L(Kj?;r)  =  -ln(detK,)-r^(K,)-1r.  (20) 

As  shown  first  by  Burg,  et  al.  [5]  a  necessary  condition  for  a  matrix  to  maximize  the  loglikelihood  function 
is  the  trace  condition 

rr[(K*>rr'K*l-K«')5K,]  =  0.  (21) 

The  matrix  8K r  is  a  variational  matrix  which  takes  values  in  all  possible  additive  variations  of  the  matrix 
K* .  The  set  of  possible  K«  is  described  by  Equation  (19).  If  we  rewrite  the  trace  condition  in  terms  of 
Kg  ,  then  we  get 

tr  [Sr  (S/Kg  St  +  NoD'H"'  -  Sf  Kg  Sr  -  NoI)(S/Kc  Sr  +  Noir'Sf  5KC  ]  =  0.  (22) 

Equation  (22)  shows  the  three  ways  in  which  the  present  spectrum  estimation  problem  differs  from  that  in 
[lj.  First,  there  is  the  signal  matrix  Sr  appearing  in  (22)  multiplying  Kg  wherever  it  appears.  Second,  the 
additive  noise  manifests  its  influence  in  the  equation  above  through  the  appearance  of  No I.  Third,  the 
matrix  Kg  in  Equation  (22)  is  a  block  diagonal  matrix  with  Toepiitz  blocks,  not  merely  Toeplitz.  Despite 
these  differences  an  EM  algorithmic  approach  to  the  solution  can  be  derived  in  much  the  same  way  as  in 
[1]. 

For  each  delay  k,  define  the  complete  data  vector  yn(k)  to  be  the  periodic  extension  of  the  data 
vector  ja(k) 

jN(k)‘  =  [jG(ky  yx(k)']  (23) 

where  7a  (k)  is  an  (N-G)xl  vector  which  augments  jo(k)  to  obtain  a  full  period  sample  of  the  periodic 
process.  Let  Jn  be  the  NIrxI  vector  made  by  stacking  the  js(k)  to  form  one  long  vector.  The  complete 
data  loglikelihood  is 

L(KH\ys)  =  -  (ln(detKv))  -  y^(Kw)-'yw,  (24) 

where  K*  is  a  block  diagonal  matrix  with  each  of  the  h  blocks  being  a  circulant  Toeplitz  matrix. 
Premultiplying  K*  by  W  and  postmultiplying  by  Wf  yields  a  diagonal  matrix.  This  diagonal  matrix 
consists  of  the  samples  of  the  spectrum,  some  of  which  are  constrained  to  be  zero  by  the  assumption  of  a 
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target  of  finite  extent  This  constraint  results  in 

K*  =W'M£*IMq,W.  (25) 

What  this  implies  is  that  after  rotating  the  data  jn  using  the  orthogonal  matrix  W,  some  of  the  entries  of  the 
resulting  vector  are  zero.  These  elements  may  be  removed  from  consideration  by  multiplying  the  resulting 
data  by  Mot  >  giving  a  vector  of  unconrelated  samples  of  a  process  whose  covariance  matrix  is  given  by  X 

c  =  Mo»Wyw  (26) 

E[cc1  =  X.  (27) 

This  vector  c  is  the  same  as  in  Equation  (17).  Under  the  assumptions  stated,  the  matrix  K s  is  not 

invertible.  The  reflectivity,  however,  almost  surely  does  not  have  a  component  in  the  null  space  of  K*  so 
we  can  make  some  sense  of  the  complete  data  loglikelihood  (24).  A  more  correct  way  to  write  the 
complete  data  loglikelihood  is  in  terms  of  the  rotated  coordinates  c  and  its  covariance  X: 

L^-t  •  <28> 

The  EM  algorithm  is  an  iterative  algorithm  which  at  each  step  updates  the  estimate  for  the  X  by 
maximizing  the  conditional  expected  value  of  the  complete  data  loglikelihood  over  X.  The  E-step  of  the 
EM  algorithm  performs  the  expected  value  of  (28)  given  the  incomplete  data  and  the  previous  estimate  for 
X.  The  M-step  consists  of  maximizing  the  result  of  this  expectation  over  the  a(i  Jc).  Since  the  complete 
data  loglikelihood  in  rotated  coordinates  separates  into  the  sum  of  independent  samples  in  (28),  the  result 
of  taking  the  maximum  over  the  spectral  values  at  step  p+1  is  just 

<j^+»(»Jk)»E[|c(tjk)j2|Z^>ar].  (29) 

The  estimate  of  the  covariance  of  y n  at  step  p+1  is  found  by  transforming  back  to  those  coordinates 

Ktr+'Kl-n  Jfc)  =  [yN(m  >NJc)  |r,K^>].  (30) 

This  equation  makes  sense  intuitively.  It  says  that  to  find  the  maximum  likelihood  estimate  over  the 
constrained  set  of  Toeplitz  covariances,  augment  the  covariance  matrix  at  each  delay  with  the  conditional 
mean  and  mean  square  estimates  of  the  missing  lags.  At  the  convergence  point  of  the  algorithm,  the 
covariance  estimates  equal  the  conditional  mean  estimates  of  the  lag  products. 

Returning  to  Equation  (27),  the  estimate  at  step  p+1  of  X  is  given  by  the  diagonal  elements  of  the 
conditional  expectation  of  ccf  or 

Z(p+t)  =  £[ccnx^4-],  (31) 

d 

where  =  means  that  the  off  diagonal  terms  on  the  left  are  zero  and  the  diagonal  terms  on  the  left  equal  the 
diagonal  terms  on  the  right  The  computation  indicated  in  (31)  is  a  standard  problem  in  estimation  theory. 
This  equation  may  be  written  as  [4] 

£o-i) = i»)  +  E»r(r'xo>r + iVoi)-1(rr^  -  no»r  -  iVoixr^r + n0 r-t  W  (32) 

Equation  (32)  defines  the  iteration  sequence  used  by  the  computations.  The  diagonal  elements  of  X  are  the 
values  which  are  displayed  as  the  scattering  function  image  of  the  target.  Thus,  at  each  stage  an  image  is 
calculated  and  may  be  displayed.  Some  appropriate  stopping  criterion  is  used  to  terminate  the  algorithm. 
It  should  be  pointed  out  that  at  each  stage  of  the  algorithm  the  conditional  mean  estimate  of  the  reflectance 
is  also  generated.  At  step  p,  this  estimate  is 

E  [cjX^^]  =  X^>nr^)r  +  WoD-'r.  (33) 

The  magnitude  or  the  magnitude  squared  of  c  are  commonly  viewed  as  images  of  the  target  Thus  both 
types  of  radar  image  commonly  viewed  are  generated  by  our  algorithm.  We  feel  this  is  a  unique  feature  of 
our  algorithm. 
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Let  the  corresponding  sequence  of  covariance  matrices  for  the  data  r  be  denoted 
K/t^  =  r'l^r  +  Since  this  iteration  is  an  EM  algorithm,  it  has  all  of  the  properties  associated  with 
this  type  of  algorithm.  In  particular,  the  incomplete  data  loglikelihood  is  nondecreasing  in  the  sequence  of 
covariance  matrices  Kj^>. 

This  section  has  presented  a  derivation  of  the  equations  used  to  produce  target  images.  This 
approach  starts  with  a  model  which  accurately  accounts  for  the  random  nature  of  radar  reflections  and 
adopts  the  maximum  likelihood  method  of  statistics  to  estimate  delay-Doppler  high  resolution  images  of 
radar  targets.  Questions  associated  with  uniqueness  of  spectrum  estimates  and  convergence  of  the  EM 
algorithm  are  discussed  in  the  following  section. 


Convergence  issues.  This  section  addresses  some  of  the  convergence  questions  associated  with  the 
the  EM  algorithms  proposed  in  earlier  sections.  These  results  are  stated  so  that  they  apply  to  both  the  radar 
imaging  problem  studied  here  and  the  Toeplitz  estimation  problem  from  [1].  Let  the  integer  M  stand  for 
l r  lot  in  the  radar  problem.  The  matrix  T  refers  to  the  T  from  the  last  section  or  simply  to  WG  if  we  are 
discussing  the  original  Toeplitz  problem.  Also,  N o  is  zero  for  the  original  Toeplitz  problem. 

One  question  of  importance  is  the  uniqueness  of  estimates  of  the  parameters  of  interest.  This 
question  is  addressed  by  looking  at  the  Cramer-Rao  bounds  on  the  variance  of  estimates.  The  Cramer-Rao 
bounds  are  obtained  by  inverting  the  Fisher  information  matrix.  When  the  Fisher  information  matrix  is 
singular,  these  bounds  are  infinite.  It  is  shown  how  singularity  of  the  Fisher  information  matrix 
corresponds  to  nonuniqueness  of  parameter  estimates. 

Definition:  Let  y*  denote  the  k*  row  of  the  A/xG  matrix  I\  The  M xG2  matrix  F  has  k*  row  given  by 
Y*  ®Y it,  where  ®  denotes  the  kronecker  product. 

Theorem  1:  The  Fisher  information  matrix  for  estimating  I  given  data  r  is  equal  to 

F(Kg{9Kft)Ft  .  (34) 

Proof:  The  Fisher  information  matrix  is  just  the  negative  of  the  expected  value  of  the  second  derivative  of 
the  log-likelihood  function.  This  second  derivative  is  evaluated  in  the  appendix  of  [4]  and  taking  the 
expected  value  yields  the  above  expression. 

Theore  m  2:  Suppose  the  Fisher  information  matrix  in  (34)  is  singular  and  that  the  matrix  Kg  is  positive 
definite.  Then  there  does  not  exist  a  £  which  is  positive  definite  which  yields  a  unique  maximum  of  the 
log-likelihood. 

Proof:  Since  the  rank  of  the  matrix  Kgx®Kgf  equals  G2,  and  by  the  form  of  the  matrix,  the  Fisher 
information  matrix  is  singular  if  and  only  if  the  matrix  F  has  rank  less  than  M  if  and  only  if  there  exists  a 
real  vector  s  such  that  Ffs  -  0.  Such  an  s  exists  if  and  only  if  there  exists  a  real  diagonal  matrix  D 
(£)  =  di  :g  (s ))  such  that 

r(Lfo£»)r=r^r  (35) 

for  all  r  ral  a.  If  £  is  positive  definite  and  maximizes  the  log- likelihood,  then  there  exists  a  {3  such  that  for 
all  0  5  a  £  p  the  matrix  £  +  aD  is  nonnegative  definite  and  yields  the  same  covariance  matrix  and  hence 
the  same  value  for  the  log-likelihood. 

CorolLry  l:  For  the  spectrum  estimation  problem  from  [1],  there  does  not  exist  a  positive  definite  £  which 
yields  a  unique  maximum  of  the  log-likelihood  if  N  >  2G  -1. 

Proof:  The  matrix  F  constructed  above  has  rank  less  than  or  equal  to  2G-1. 

Note  that  this  theorem  does  not  say  that  the  estimate  of  the  Toeplitz  covariance  matrices  generated 
by  the  algorithm  are  not  unique.  The  theorem  and  its  corollary  relate  to  the  uniqueness  of  the  spectrum 
samples.  For  some  problems  the  parameters  of  interest  are  in  the  covariance  matrix  K*  or  in  the  Toeplitz 
matrix  Kg.  There  could  be  (and  indeed  are  for  M  large  enough)  many  £  which  yield  the  same  estimate  for 
Kg .  Theorem  2  can  be  applied  to  problems  where  one  desires  to  know  how  big  to  make  N .  In  general,  it 
is  desired  to  have  N  as  short  as  possible  in  order  to  reduce  the  number  of  parameters  to  be  estimated.  If  it 
can  be  shown  that  with  a  given  N  there  exists  a  positive  definite  matrix  which  yields  the  maximum 
likelihood  estimate  for  the  covariance  matrix  and  the  conditions  of  the  theorem  are  satisfied,  then  one 
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might  consider  using  a  smaller  N  to  reduce  the  size  of  Z. 

For  some  problems,  including  the  radar  problem,  it  is  the  matrix  I  which  is  of  interest.  For  these 
problems,  it  is  very  important  to  know  when  the  estimate  is  unique.  The  radar  imaging  problem  is  one 
example.  In  the  radar  problem,  the  result  which  is  displayed  as  an  image  is  an  array  whose  elements  are 
the  diagonal  entries  from  £.  In  order  to  be  able  to  generate  a  unique  image,  the  conditions  of  theorem  2 
must  be  satisfied. 

Some  of  the  issues  associated  with  convergence  of  the  EM  algorithm  for  the  problems  described  are 
addressed  next  The  following  material  is  adapted  from  the  material  in  [1]  to  be  applicable  to  the  radar 
problem  as  well. 

Definition:  Let  Kg  be  the  set  of  positive  definite  matrices  K*  given  by  Equation  (19)  whose  entries  are 
bounded  by  some  number  a .  Let  Kg  be  the  closure  of  Kg ,  the  set  of  nonnegative  definite  matrices  of  this 
parameterized  form. 

The  important  issues  for  the  following  theorems  are  not  the  matrices  which  go  into  F.  What  is 
important  is  that  any  Kg  e  Kg  may  be  written  as 


where  y*  is  the  k*  row  of  T  which  is  fixed  once  the  model  is  specified.  The  only  parameters  which  must 
be  found  are  the  a(k)  which  are  specified  to  be  greater  than  or  equal  to  zero.  An  element  of  Kg  must  have 
,V0  equal  to  zero.  This  is  an  important  cas£  for  which  we  wish  to  guarantee  that  the  estimate  of  K«  is 
nonsingular.  Clearly  for  any  fixed  b  the  set  Kg  is  compact 

Theorem  3:  Let  T  be  any  fixed  M  xG  matrix  with  complex  entries.  Let  r  be  an  observation  of  a  Gy.  1  0- 
mean  Gaussian  random  vector  whose  covariance  is  some  positive  definite  hermitian  symmetric  matrix. 
Then 

a)  There  does  not  exist  a  singular  Kg  e  Kg  such  that  r  is  in  the  range  space  of  Kg ,  with  probability  one. 

b)  The  log-likelihood  function  is  bounded  from  above  over  the  set  Kg,  with  probability  one. 

Proof:  a)  Suppose  that  any  G  rows  of  T  are  linearly  independent  Since  K*  is  given  by  (36),  a  singular 
matrix  in  this  class  must  be  given  by 

Ki=£a(k)y£yk,  (37) 

where  Jc{ 0,l,2,...,Af-l/  consists  of  G-l  or  fewer  integers  which  correspond  to  the  nonzero  diagonal 
entries  of  E.  Since  the  true  covariance  for  r  is  nonsingular,  the  probability  that  r  lies  in  the  subspace 
spanned  by  (y{\ keJ}  is  zero.  Since  there  are  a  finite  number  of  such  spaces  and  the  probability  that  r  is  in 
any  one  of  them  is  zero,  the  probability  that  r  is  in  the  range  of  any  singular  K*  in  the  set  is  zero.  If  any  G 
rows  are  not  independent,  then  singular  matrices  may  be  written  as  the  sum  of  more  than  G-l  outer 
products  ytjk-  But  the  data  would  still  have  to  lie  in  a  subspace  spanned  by  fewer  than  G  independent 
vectors  and  thus  the  probability  of  this  is  zero  and  this  part  of  the  theorem  follows, 
b)  This  part  follows  from  [6]  where  it  is  shown  that  the  log-likelihood  function  is  bounded  above  when  the 
Harq  are  not  in  the  range  space  of  a  singular  covariance  matrix  in  the  set  in  question.  The  proof  is  based  on 
the  following  facts.  First,  if  Kg  is  nonsingular  and  its  eigenvalues  are  bounded  from  above  and  below,  the 
log-likelihood  is  bounded.  Second,  if  K#  is  singular  and  the  data  are  in  its  range,  the  log-likelihood  is 
unbounded  above;  but  this  is  a  zero  probability  event.  Third,  if  K*  is  singular  (with  rank  n)  and  the  data  are 
not  in  its  range,  the  log-likelihood  is  unbounded  from  below.  This  is  shown  by  writing  Kg  as  the  limit  as 
e  -►  0  of  K#  +  F(el)r.  We  examine  the  loglikelihood  (20)  when  r  is  not  in  the  range  of  Kg .  The  term 
-ln(det(K#  +  eT^T))  has  bounded  terms  and  an  unbounded  component  of  the  form  -(G  -n  )ln(e).  The  term 
-r^(K*  +  eT*T)-lr  has  bounded  terms  and  an  unbounded  component  of  the  form  -c/e,  where  c  is  some 
positive  number.  Then, 

linj  — (G  — n  )ln(c) —  =  ~°°  •  (38) 

There  are  more  facts  needed  to  prove  convergence.  One  fact  is  that  the  log-likelihood  is  increasing 
at  each  step  of  the  algorithm.  Another  is  that  the  iterates  stay  in  a  bounded  set  so  that  the  above  theorems 
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apply.  The  theorems  before  this  point  apply  to  the  problem  no  matter  what  algorithm  is  used  to  find  the 
maximum  likelihood  estimate,  while  the  theorems  that  follow  apply  for  our  particular  algorithm. 

Theorem  4:  The  iterates  defined  by  the  EM  algorithm  (32)  produce  a  sequence  of  log- likelihoods  which 
are  nondecreasing, 

L (Kjf +1>;r)  - L (K^>;r)  > Q (K^+‘> I K^>) -  Q  W I W  5  0  , 
where  L(.;.)is  the  log-likelihood  for  the  problem. 

Proof:  This  is  just  a  result  of  the  sequence  being  generated  by  an  EM  algorithm  [7, 8]. 

Theorem  5:  L  (K^y)  =  L  (K£»;r)  if  and  only  =  Z»>. 

Proof:  This  is  a  result  of  the  concavity  of  the  complete  data  log-likelihood.  Take  the  second  derivative 
with  respect  to  the  variable  being  maximized,  Z.  For  each  diagonal  entry  of  Z  this  derivative  is  either 
positive  or  zero.  It’s  zero  if  and  only  if  the  previous  corresponding  entry  of  Z  is  zero,  and  this  entry  would 
remain  zero.  Thus  the  maximizer  is  unique  and  is  given  by  Z^*l).  By  the  inequality  from  the  theorem  3, 
theorem  4  follows. 

The  one  last  theorem  we  would  like  to  have  is  that  the  iterates  remain  in  a  bounded  set.  It  has  been 
our  experience  in  computations  that  the  iterates  do  remain  bounded,  but  we  have  had  trouble  proving  this  in 
the  general  case.  We  have  observed  in  computations  that  Z  may  tend  to  a  singular  limit  This  is  not 
precluded  by  any  of  the  above  theorems.  In  fact,  for  our  radar  imaging  problem  we  do  not  wish  to  exclude 
this  possibility  since  a  zero  estimate  of  the  power  reflected  from  a  point  simply  means  that  there  is  no  target 
at  that  point 


Conclusions.  We  have  presented  an  algorithm  for  generating  images  of  radar  targets  in  the  delay- 
Doppler  plane.  The  approach  has  been  estimation  based  because  of  our  assumption  of  targets  which  are 
rough  on  the  order  of  a  wavelength  of  the  carrier  frequency.  Some  of  the  theoretical  properties  of  this 
approach  and  uniqueness  of  estimates  have  been  discussed.  Presently  we  are  implementing  the  proposed 
algorithm  and  performing  computational  studies.  We  are  also  addressing  several  theoretical  issues.  One 
issue  of  particular  importance  is  the  incorporation  of  specular  components  in  the  algorithm.  These  points 
would  have  a  different  statistical  characterization  than  the  diffuse  components  considered  here  and  a 
correspondingly  altered  loglikelihood  to  be  maximized.  Computationally,  we  are  examining  the 
convergence  of  our  algorithm,  its  computational  complexity,  and  comparing  the  performance  of  the 
algorithm  to  other  approaches. 
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